MOM6
Data Types | Functions/Subroutines
midas_vertmap Module Reference

Detailed Description

Routines for initialization callable from MOM6 or Python (MIDAS)

Data Types

interface  fill_boundaries
 Fill grid edges. More...
 

Functions/Subroutines

real function, dimension(size(tr_in, 1), size(tr_in, 2), nlay), public tracer_z_init (tr_in, z_edges, e, nkml, nkbl, land_fill, wet, nlay, nlevs, debug, i_debug, j_debug, eps_z)
 Layer model routine for remapping tracers. More...
 
integer function, dimension(size(a, 1), size(x, 1)) bisect_fast (a, x, lo, hi)
 Return the index where to insert item x in list a, assuming a is sorted. The return values [i] is such that all e in a[:i-1] have e <= x, and all e in a[i:] have e > x. So if x already appears in the list, will insert just after the rightmost x already there. Optional args lo (default 1) and hi (default len(a)) bound the slice of a to be searched. More...
 
subroutine, public determine_temperature (temp, salt, R, p_ref, niter, land_fill, h, k_start, eos)
 This subroutine determines the potential temperature and salinity that is consistent with the target density using provided initial guess. More...
 
subroutine find_overlap (e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z2)
 This subroutine determines the layers bounded by interfaces e that overlap with the depth range between Z_top and Z_bot, and also the fractional weights of each layer. It also calculates the normalized relative depths of the range of each layer that overlaps that depth range. Note that by convention, e decreases with increasing k and Z_top > Z_bot. More...
 
real function find_limited_slope (val, e, k)
 This subroutine determines a limited slope for val to be advected with a piecewise limited scheme. More...
 
real function, dimension(size(rho, 1), size(rho, 2), size(rb, 1)), public find_interfaces (rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps_z, eps_rho)
 Find interface positions corresponding to density profile. More...
 
subroutine, public meshgrid (x, y, x_T, y_T)
 Create a 2d-mesh of grid coordinates from 1-d arrays. More...
 
subroutine smooth_heights (zi, fill, bad, sor, niter, cyclic_x, tripolar_n)
 Solve del2 (zi) = 0 using successive iterations with a 5 point stencil. Only points fill==1 are modified. Except where bad==1, information propagates isotropically in index space. The resulting solution in each region is an approximation to del2(zi)=0 subject to boundary conditions along the valid points curve bounding this region. More...
 
integer function, dimension(0:size(m, 1)+1, 0:size(m, 2)+1) fill_boundaries_int (m, cyclic_x, tripolar_n)
 Fill grid edges. More...
 
real function, dimension(0:size(m, 1)+1, 0:size(m, 2)+1) fill_boundaries_real (m, cyclic_x, tripolar_n)
 fill grid edges More...
 

Function/Subroutine Documentation

◆ bisect_fast()

integer function, dimension(size(a,1),size(x,1)) midas_vertmap::bisect_fast ( real, dimension(:,:), intent(in)  a,
real, dimension(:), intent(in)  x,
integer, dimension(size(a,1)), intent(in), optional  lo,
integer, dimension(size(a,1)), intent(in), optional  hi 
)
private

Return the index where to insert item x in list a, assuming a is sorted. The return values [i] is such that all e in a[:i-1] have e <= x, and all e in a[i:] have e > x. So if x already appears in the list, will insert just after the rightmost x already there. Optional args lo (default 1) and hi (default len(a)) bound the slice of a to be searched.

Parameters
[in]aSorted list
[in]xItem to be inserted
[in]loLower bracket of optional range to search
[in]hiUpper bracket of optional range to search

Definition at line 302 of file midas_vertmap.F90.

302  real, dimension(:,:), intent(in) :: a !< Sorted list
303  real, dimension(:), intent(in) :: x !< Item to be inserted
304  integer, dimension(size(a,1)), optional, intent(in) :: lo !< Lower bracket of optional range to search
305  integer, dimension(size(a,1)), optional, intent(in) :: hi !< Upper bracket of optional range to search
306  integer, dimension(size(a,1),size(x,1)) :: bi_r
307 
308  integer :: mid,num_x,num_a,i
309  integer, dimension(size(a,1)) :: lo_,hi_,lo0,hi0
310  integer :: nprofs,j
311 
312  lo_=1;hi_=size(a,2);num_x=size(x,1);bi_r=-1;nprofs=size(a,1)
313 
314  if (PRESENT(lo)) then
315  where (lo>0) lo_=lo
316  endif
317  if (PRESENT(hi)) then
318  where (hi>0) hi_=hi
319  endif
320 
321  lo0=lo_;hi0=hi_
322 
323  do j=1,nprofs
324  do i=1,num_x
325  lo_=lo0;hi_=hi0
326  do while (lo_(j) < hi_(j))
327  mid = (lo_(j)+hi_(j))/2
328  if (x(i) < a(j,mid)) then
329  hi_(j) = mid
330  else
331  lo_(j) = mid+1
332  endif
333  enddo
334  bi_r(j,i)=lo_(j)
335  enddo
336  enddo
337 
338 
339  return
340 

Referenced by find_interfaces().

Here is the caller graph for this function:

◆ determine_temperature()

subroutine, public midas_vertmap::determine_temperature ( real, dimension(:,:,:), intent(inout)  temp,
real, dimension(:,:,:), intent(inout)  salt,
real, dimension(size(temp,3)), intent(in)  R,
real, intent(in)  p_ref,
integer, intent(in)  niter,
real, intent(in)  land_fill,
real, dimension(:,:,:), intent(in)  h,
integer, intent(in)  k_start,
type(eos_type), pointer  eos 
)

This subroutine determines the potential temperature and salinity that is consistent with the target density using provided initial guess.

Parameters
[in,out]temppotential temperature [degC]
[in,out]saltsalinity [PSU]
[in]rdesired potential density [kg m-3].
[in]p_refreference pressure [Pa].
[in]nitermaximum number of iterations
[in]k_startstarting index (i.e. below the buffer layer)
[in]land_fillland fill value
[in]hlayer thickness, used only to avoid working on massless layers
eosseawater equation of state control structure

Definition at line 364 of file midas_vertmap.F90.

364  real, dimension(:,:,:), intent(inout) :: temp !< potential temperature [degC]
365  real, dimension(:,:,:), intent(inout) :: salt !< salinity [PSU]
366  real, dimension(size(temp,3)), intent(in) :: R !< desired potential density [kg m-3].
367  real, intent(in) :: p_ref !< reference pressure [Pa].
368  integer, intent(in) :: niter !< maximum number of iterations
369  integer, intent(in) :: k_start !< starting index (i.e. below the buffer layer)
370  real, intent(in) :: land_fill !< land fill value
371  real, dimension(:,:,:), intent(in) :: h !< layer thickness, used only to avoid working on massless layers
372  type(eos_type), pointer :: eos !< seawater equation of state control structure
373 
374  real, parameter :: T_max = 31.0, t_min = -2.0
375 #endif
376  ! Local variables (All of which need documentation!)
377  real(kind=8), dimension(size(temp,1),size(temp,3)) :: t, s, dt, ds, rho, hin
378  real(kind=8), dimension(size(temp,1),size(temp,3)) :: drho_dt, drho_ds
379  real(kind=8), dimension(size(temp,1)) :: press
380  integer :: nx, ny, nz, nt, i, j, k, n, itt
381  real :: dT_dS
382  logical :: adjust_salt, old_fit
383  real, parameter :: S_min = 0.5, s_max=65.0
384  real, parameter :: tol=1.e-4, max_t_adj=1.0, max_s_adj = 0.5
385 
386  old_fit = .true. ! reproduces siena behavior
387  ! will switch to the newer method which simultaneously adjusts
388  ! temp and salt based on the ratio of the thermal and haline coefficients.
389 
390  nx=size(temp,1) ; ny=size(temp,2) ; nz=size(temp,3)
391 
392  press(:) = p_ref
393 
394  do j=1,ny
395  ds(:,:) = 0. ! Needs to be zero everywhere since there is a maxval(abs(dS)) later...
396  t=temp(:,j,:)
397  s=salt(:,j,:)
398  hin=h(:,j,:)
399  dt=0.0
400  adjust_salt = .true.
401  iter_loop: do itt = 1,niter
402 #ifdef PY_SOLO
403  rho=wright_eos_2d(t,s,p_ref)
404  drho_dt=alpha_wright_eos_2d(t,s,p_ref)
405 #else
406  do k=1, nz
407  call calculate_density(t(:,k), s(:,k), press, rho(:,k), 1, nx, eos)
408  call calculate_density_derivs(t(:,k), s(:,k), press, drho_dt(:,k), drho_ds(:,k), 1, nx, eos)
409  enddo
410 #endif
411  do k=k_start,nz ; do i=1,nx
412 
413 ! if (abs(rho(i,k)-R(k))>tol .and. hin(i,k)>epsln .and. abs(T(i,k)-land_fill) < epsln) then
414  if (abs(rho(i,k)-r(k))>tol) then
415  if (old_fit) then
416  dt(i,k) = max(min((r(k)-rho(i,k)) / drho_dt(i,k), max_t_adj), -max_t_adj)
417  t(i,k) = max(min(t(i,k)+dt(i,k), t_max), t_min)
418  else
419  dt_ds = 10.0 - min(-drho_dt(i,k)/drho_ds(i,k),10.)
420  !### RWH: Based on the dimensions alone, the expression above should be:
421  ! dT_dS = 10.0 - min(-drho_dS(i,k)/drho_dT(i,k),10.)
422  ds(i,k) = (r(k)-rho(i,k)) / (drho_ds(i,k) - drho_dt(i,k)*dt_ds )
423  dt(i,k) = -dt_ds*ds(i,k)
424  ! dT(i,k) = max(min(dT(i,k), max_t_adj), -max_t_adj)
425  t(i,k) = max(min(t(i,k)+dt(i,k), t_max), t_min)
426  s(i,k) = max(min(s(i,k)+ds(i,k), s_max), s_min)
427  endif
428  endif
429  enddo ; enddo
430  if (maxval(abs(dt)) < tol) then
431  adjust_salt = .false.
432  exit iter_loop
433  endif
434  enddo iter_loop
435 
436  if (adjust_salt .and. old_fit) then ; do itt = 1,niter
437 #ifdef PY_SOLO
438  rho = wright_eos_2d(t,s,p_ref)
439  drho_ds = beta_wright_eos_2d(t,s,p_ref)
440 #else
441  do k=1, nz
442  call calculate_density(t(:,k),s(:,k),press,rho(:,k),1,nx,eos)
443  call calculate_density_derivs(t(:,k),s(:,k),press,drho_dt(:,k),drho_ds(:,k),1,nx,eos)
444  enddo
445 #endif
446  do k=k_start,nz ; do i=1,nx
447 ! if (abs(rho(i,k)-R(k))>tol .and. hin(i,k)>epsln .and. abs(T(i,k)-land_fill) < epsln ) then
448  if (abs(rho(i,k)-r(k)) > tol) then
449  ds(i,k) = max(min((r(k)-rho(i,k)) / drho_ds(i,k), max_s_adj), -max_s_adj)
450  s(i,k) = max(min(s(i,k)+ds(i,k), s_max), s_min)
451  endif
452  enddo ; enddo
453  if (maxval(abs(ds)) < tol) exit
454  enddo ; endif
455 
456  temp(:,j,:)=t(:,:)
457  salt(:,j,:)=s(:,:)
458  enddo
459 

◆ fill_boundaries_int()

integer function, dimension(0:size(m,1)+1,0:size(m,2)+1) midas_vertmap::fill_boundaries_int ( integer, dimension(:,:), intent(in)  m,
logical, intent(in)  cyclic_x,
logical, intent(in)  tripolar_n 
)
private

Fill grid edges.

Parameters
[in]minput array
[in]cyclic_xzonal cyclic condition
[in]tripolar_nnorthern fold condition
Returns
output filled array

Definition at line 782 of file midas_vertmap.F90.

782  integer, dimension(:,:), intent(in) :: m !< input array
783  logical, intent(in) :: cyclic_x !< zonal cyclic condition
784  logical, intent(in) :: tripolar_n !< northern fold condition
785  integer, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp !< output filled array
786  ! Local variables
787  real, dimension(size(m,1),size(m,2)) :: m_real
788  real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real
789 
790  m_real = real(m)
791 
792  mp_real = fill_boundaries_real(m_real,cyclic_x,tripolar_n)
793 
794  mp = int(mp_real)
795 
796  return
797 

References fill_boundaries_real().

Here is the call graph for this function:

◆ fill_boundaries_real()

real function, dimension(0:size(m,1)+1,0:size(m,2)+1) midas_vertmap::fill_boundaries_real ( real, dimension(:,:), intent(in)  m,
logical, intent(in)  cyclic_x,
logical, intent(in)  tripolar_n 
)
private

fill grid edges

Parameters
[in]minput array
[in]cyclic_xzonal cyclic condition
[in]tripolar_nnorthern fold condition
Returns
output filled array

Definition at line 802 of file midas_vertmap.F90.

802  real, dimension(:,:), intent(in) :: m !< input array
803  logical, intent(in) :: cyclic_x !< zonal cyclic condition
804  logical, intent(in) :: tripolar_n !< northern fold condition
805  real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp !< output filled array
806 
807  integer :: ni,nj,i,j
808 
809  ni=size(m,1); nj=size(m,2)
810 
811  mp(1:ni,1:nj)=m(:,:)
812 
813  if (cyclic_x) then
814  mp(0,1:nj)=m(ni,1:nj)
815  mp(ni+1,1:nj)=m(1,1:nj)
816  else
817  mp(0,1:nj)=m(1,1:nj)
818  mp(ni+1,1:nj)=m(ni,1:nj)
819  endif
820 
821  mp(1:ni,0)=m(1:ni,1)
822  if (tripolar_n) then
823  do i=1,ni
824  mp(i,nj+1)=m(ni-i+1,nj)
825  enddo
826  else
827  mp(1:ni,nj+1)=m(1:ni,nj)
828  endif
829 
830  return
831 

Referenced by fill_boundaries_int().

Here is the caller graph for this function:

◆ find_interfaces()

real function, dimension(size(rho,1),size(rho,2),size(rb,1)), public midas_vertmap::find_interfaces ( real, dimension(:,:,:), intent(in)  rho,
real, dimension(size(rho,3)), intent(in)  zin,
real, dimension(:), intent(in)  Rb,
real, dimension(size(rho,1),size(rho,2)), intent(in)  depth,
real, dimension(size(rho,1),size(rho,2)), intent(in), optional  nlevs,
integer, intent(in), optional  nkml,
integer, intent(in), optional  nkbl,
real, intent(in), optional  hml,
logical, intent(in), optional  debug,
real, intent(in), optional  eps_z,
real, intent(in), optional  eps_rho 
)

Find interface positions corresponding to density profile.

Parameters
[in]rhopotential density in z-space [kg m-3 or R ~> kg m-3]
[in]zinInput data levels [m or Z ~> m].
[in]rbtarget interface densities [kg m-3 or R ~> kg m-3]
[in]depthocean depth [Z ~> m].
[in]nlevsnumber of valid points in each column
[in]debugoptional debug flag
[in]nkmlnumber of mixed layer pieces
[in]nkblnumber of buffer layer pieces
[in]hmlmixed layer depth [Z ~> m].
[in]eps_zA negligibly small layer thickness [m or Z ~> m].
[in]eps_rhoA negligibly small density difference [kg m-3 or R ~> kg m-3].
Returns
The returned interface, in the same units az zin.

Definition at line 563 of file midas_vertmap.F90.

563  real, dimension(:,:,:), &
564  intent(in) :: rho !< potential density in z-space [kg m-3 or R ~> kg m-3]
565  real, dimension(size(rho,3)), &
566  intent(in) :: zin !< Input data levels [m or Z ~> m].
567  real, dimension(:), intent(in) :: Rb !< target interface densities [kg m-3 or R ~> kg m-3]
568  real, dimension(size(rho,1),size(rho,2)), &
569  intent(in) :: depth !< ocean depth [Z ~> m].
570  real, dimension(size(rho,1),size(rho,2)), &
571  optional, intent(in) :: nlevs !< number of valid points in each column
572  logical, optional, intent(in) :: debug !< optional debug flag
573  integer, optional, intent(in) :: nkml !< number of mixed layer pieces
574  integer, optional, intent(in) :: nkbl !< number of buffer layer pieces
575  real, optional, intent(in) :: hml !< mixed layer depth [Z ~> m].
576  real, optional, intent(in) :: eps_z !< A negligibly small layer thickness [m or Z ~> m].
577  real, optional, intent(in) :: eps_rho !< A negligibly small density difference [kg m-3 or R ~> kg m-3].
578  real, dimension(size(rho,1),size(rho,2),size(Rb,1)) :: zi !< The returned interface, in the same units az zin.
579 
580  ! Local variables
581  real, dimension(size(rho,1),size(rho,3)) :: rho_ ! A slice of densities [R ~> kg m-3]
582  real, dimension(size(rho,1)) :: depth_
583  logical :: unstable
584  integer :: dir
585  integer, dimension(size(rho,1),size(Rb,1)) :: ki_
586  real, dimension(size(rho,1),size(Rb,1)) :: zi_
587  integer, dimension(size(rho,1),size(rho,2)) :: nlevs_data
588  integer, dimension(size(rho,1)) :: lo, hi
589  real :: slope,rsm,drhodz,hml_
590  integer :: n,i,j,k,l,nx,ny,nz,nt
591  integer :: nlay,kk,nkml_,nkbl_
592  logical :: debug_ = .false.
593  real :: epsln_Z ! A negligibly thin layer thickness [m or Z ~> m].
594  real :: epsln_rho ! A negligibly small density change [kg m-3 or R ~> kg m-3].
595  real, parameter :: zoff=0.999
596 
597  nlay=size(rb)-1
598 
599  zi(:,:,:) = 0.0
600 
601  if (PRESENT(debug)) debug_=debug
602 
603  nx = size(rho,1); ny=size(rho,2); nz = size(rho,3)
604  nlevs_data(:,:) = size(rho,3)
605 
606  nkml_ = 0 ; if (PRESENT(nkml)) nkml_ = max(0, nkml)
607  nkbl_ = 0 ; if (PRESENT(nkbl)) nkbl_ = max(0, nkbl)
608  hml_ = 0.0 ; if (PRESENT(hml)) hml_ = hml
609  epsln_z = 1.0e-10 ; if (PRESENT(eps_z)) epsln_z = eps_z
610  epsln_rho = 1.0e-10 ; if (PRESENT(eps_rho)) epsln_rho = eps_rho
611 
612  if (PRESENT(nlevs)) then
613  nlevs_data(:,:) = nlevs(:,:)
614  endif
615 
616  do j=1,ny
617  rho_(:,:) = rho(:,j,:)
618  i_loop: do i=1,nx
619  if (debug_) then
620  print *,'looking for interfaces, i,j,nlevs= ',i,j,nlevs_data(i,j)
621  print *,'initial density profile= ', rho_(i,:)
622  endif
623  unstable=.true.
624  dir=1
625  do while (unstable)
626  unstable=.false.
627  if (dir == 1) then
628  do k=2,nlevs_data(i,j)-1
629  if (rho_(i,k) - rho_(i,k-1) < 0.0 ) then
630  if (k == 2) then
631  rho_(i,k-1) = rho_(i,k)-epsln_rho
632  else
633  drhodz = (rho_(i,k+1)-rho_(i,k-1)) / (zin(k+1)-zin(k-1))
634  if (drhodz < 0.0) unstable=.true.
635  rho_(i,k) = rho_(i,k-1) + drhodz*zoff*(zin(k)-zin(k-1))
636  endif
637  endif
638  enddo
639  dir = -1*dir
640  else
641  do k=nlevs_data(i,j)-1,2,-1
642  if (rho_(i,k+1) - rho_(i,k) < 0.0) then
643  if (k == nlevs_data(i,j)-1) then
644  rho_(i,k+1) = rho_(i,k-1)+epsln_rho
645  else
646  drhodz = (rho_(i,k+1)-rho_(i,k-1))/(zin(k+1)-zin(k-1))
647  if (drhodz < 0.0) unstable=.true.
648  rho_(i,k) = rho_(i,k+1)-drhodz*(zin(k+1)-zin(k))
649  endif
650  endif
651  enddo
652  dir = -1*dir
653  endif
654  enddo
655  if (debug_) then
656  print *,'final density profile= ', rho_(i,:)
657  endif
658  enddo i_loop
659 
660  ki_(:,:) = 0
661  zi_(:,:) = 0.0
662  depth_(:) = -1.0*depth(:,j)
663  lo(:) = 1
664  hi(:) = nlevs_data(:,j)
665  ki_ = bisect_fast(rho_, rb, lo, hi)
666  ki_(:,:) = max(1, ki_(:,:)-1)
667  do i=1,nx
668  do l=2,nlay
669  slope = (zin(ki_(i,l)+1) - zin(ki_(i,l))) / max(rho_(i,ki_(i,l)+1) - rho_(i,ki_(i,l)),epsln_rho)
670  zi_(i,l) = -1.0*(zin(ki_(i,l)) + slope*(rb(l)-rho_(i,ki_(i,l))))
671  zi_(i,l) = max(zi_(i,l), depth_(i))
672  zi_(i,l) = min(zi_(i,l), -1.0*hml_)
673  enddo
674  zi_(i,nlay+1) = depth_(i)
675  do l=2,nkml_+1
676  zi_(i,l) = max(hml_*((1.0-real(l))/real(nkml_)), depth_(i))
677  enddo
678  do l=nlay,nkml_+2,-1
679  if (zi_(i,l) < zi_(i,l+1) + epsln_z) zi_(i,l) = zi_(i,l+1) + epsln_z
680  if (zi_(i,l) > -1.0*hml_) zi_(i,l) = max(-1.0*hml_, depth_(i))
681  enddo
682  enddo
683  zi(:,j,:) = zi_(:,:)
684  enddo
685 

References bisect_fast().

Here is the call graph for this function:

◆ find_limited_slope()

real function midas_vertmap::find_limited_slope ( real, dimension(:), intent(in)  val,
real, dimension(:), intent(in)  e,
integer, intent(in)  k 
)
private

This subroutine determines a limited slope for val to be advected with a piecewise limited scheme.

Parameters
[in]valAn column the values that are being interpolated.
[in]eA column's interface heights [Z ~> m] or other units.
[in]kThe layer whose slope is being determined.
Returns
The normalized slope in the intracell distribution of val.

Definition at line 530 of file midas_vertmap.F90.

530  real, dimension(:), intent(in) :: val !< An column the values that are being interpolated.
531  real, dimension(:), intent(in) :: e !< A column's interface heights [Z ~> m] or other units.
532  integer, intent(in) :: k !< The layer whose slope is being determined.
533  real :: slope !< The normalized slope in the intracell distribution of val.
534  ! Local variables
535  real :: amn, cmn
536  real :: d1, d2
537 
538  if ((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) then
539  slope = 0.0 ! ; curvature = 0.0
540  else
541  d1 = 0.5*(e(k-1)-e(k+1)) ; d2 = 0.5*(e(k)-e(k+2))
542  if (d1*d2 > 0.0) then
543  slope = ((d1**2)*(val(k+1) - val(k)) + (d2**2)*(val(k) - val(k-1))) * &
544  (e(k) - e(k+1)) / (d1*d2*(d1+d2))
545  ! slope = 0.5*(val(k+1) - val(k-1))
546  ! This is S.J. Lin's form of the PLM limiter.
547  amn = min(abs(slope), 2.0*(max(val(k-1), val(k), val(k+1)) - val(k)))
548  cmn = 2.0*(val(k) - min(val(k-1), val(k), val(k+1)))
549  slope = sign(1.0, slope) * min(amn, cmn)
550 
551  ! min(abs(slope), 2.0*(max(val(k-1),val(k),val(k+1)) - val(k)), &
552  ! 2.0*(val(k) - min(val(k-1),val(k),val(k+1))))
553  ! curvature = 0.0
554  else
555  slope = 0.0 ! ; curvature = 0.0
556  endif
557  endif
558 

Referenced by tracer_z_init().

Here is the caller graph for this function:

◆ find_overlap()

subroutine midas_vertmap::find_overlap ( real, dimension(:), intent(in)  e,
real, intent(in)  Z_top,
real, intent(in)  Z_bot,
integer, intent(in)  k_max,
integer, intent(in)  k_start,
integer, intent(out)  k_top,
integer, intent(out)  k_bot,
real, dimension(:), intent(out)  wt,
real, dimension(:), intent(out)  z1,
real, dimension(:), intent(out)  z2 
)
private

This subroutine determines the layers bounded by interfaces e that overlap with the depth range between Z_top and Z_bot, and also the fractional weights of each layer. It also calculates the normalized relative depths of the range of each layer that overlaps that depth range. Note that by convention, e decreases with increasing k and Z_top > Z_bot.

Parameters
[in]eThe interface positions, [Z ~> m] or other units.
[in]z_topThe top of the range being mapped to, [Z ~> m] or other units.
[in]z_botThe bottom of the range being mapped to, [Z ~> m] or other units.
[in]k_maxThe number of valid layers.
[in]k_startThe layer at which to start searching.
[out]k_topThe index of the top layer that overlap with the depth range.
[out]k_botThe index of the bottom layer that overlap with the depth range.
[out]wtThe relative weights of each layer from k_top to k_bot [nondim].
[out]z1Depth of the top limit of layer that contributes to a level [nondim].
[out]z2Depth of the bottom limit of layer that contributes to a level [nondim].

Definition at line 468 of file midas_vertmap.F90.

468  real, dimension(:), intent(in) :: e !< The interface positions, [Z ~> m] or other units.
469  real, intent(in) :: Z_top !< The top of the range being mapped to, [Z ~> m] or other units.
470  real, intent(in) :: Z_bot !< The bottom of the range being mapped to, [Z ~> m] or other units.
471  integer, intent(in) :: k_max !< The number of valid layers.
472  integer, intent(in) :: k_start !< The layer at which to start searching.
473  integer, intent(out) :: k_top !< The index of the top layer that overlap with the depth range.
474  integer, intent(out) :: k_bot !< The index of the bottom layer that overlap with the depth range.
475  real, dimension(:), intent(out) :: wt !< The relative weights of each layer from k_top to k_bot [nondim].
476  real, dimension(:), intent(out) :: z1 !< Depth of the top limit of layer that contributes to a level [nondim].
477  real, dimension(:), intent(out) :: z2 !< Depth of the bottom limit of layer that contributes to a level [nondim].
478 
479  ! Local variables
480  real :: Ih, e_c, tot_wt, I_totwt
481  integer :: k
482 
483  wt(:)=0.0 ; z1(:)=0.0 ; z2(:)=0.0
484  k_top = k_start ; k_bot = k_start ; wt(1) = 1.0 ; z1(1) = -0.5 ; z2(1) = 0.5
485 
486  do k=k_start,k_max ; if (e(k+1) < z_top) exit ; enddo
487  k_top = k
488 
489  if (k>k_max) return
490 
491  ! Determine the fractional weights of each layer.
492  ! Note that by convention, e and Z_int decrease with increasing k.
493  if (e(k+1) <= z_bot) then
494  wt(k) = 1.0 ; k_bot = k
495  ih = 0.0 ; if (e(k) /= e(k+1)) ih = 1.0 / (e(k)-e(k+1))
496  e_c = 0.5*(e(k)+e(k+1))
497  z1(k) = (e_c - min(e(k), z_top)) * ih
498  z2(k) = (e_c - z_bot) * ih
499  else
500  wt(k) = min(e(k),z_top) - e(k+1) ; tot_wt = wt(k) ! These are always > 0.
501  ! Ih = 0.0 ; if (e(K) /= e(K+1)) Ih = 1.0 / (e(K)-e(K+1))
502  if (e(k) /= e(k+1)) then
503  z1(k) = (0.5*(e(k)+e(k+1)) - min(e(k), z_top)) / (e(k)-e(k+1))
504  else ; z1(k) = -0.5 ; endif
505  z2(k) = 0.5
506  k_bot = k_max
507  do k=k_top+1,k_max
508  if (e(k+1) <= z_bot) then
509  k_bot = k
510  wt(k) = e(k) - z_bot ; z1(k) = -0.5
511  if (e(k) /= e(k+1)) then
512  z2(k) = (0.5*(e(k)+e(k+1)) - z_bot) / (e(k)-e(k+1))
513  else ; z2(k) = 0.5 ; endif
514  else
515  wt(k) = e(k) - e(k+1) ; z1(k) = -0.5 ; z2(k) = 0.5
516  endif
517  tot_wt = tot_wt + wt(k) ! wt(k) is always > 0.
518  if (k>=k_bot) exit
519  enddo
520 
521  i_totwt = 0.0 ; if (tot_wt > 0.0) i_totwt = 1.0 / tot_wt
522  do k=k_top,k_bot ; wt(k) = i_totwt*wt(k) ; enddo
523  endif
524 

Referenced by tracer_z_init().

Here is the caller graph for this function:

◆ meshgrid()

subroutine, public midas_vertmap::meshgrid ( real, dimension(:), intent(in)  x,
real, dimension(:), intent(in)  y,
real, dimension(size(x,1),size(y,1)), intent(inout)  x_T,
real, dimension(size(x,1),size(y,1)), intent(inout)  y_T 
)

Create a 2d-mesh of grid coordinates from 1-d arrays.

Parameters
[in]xinput x coordinates
[in]yinput y coordinates
[in,out]x_toutput 2-d version
[in,out]y_toutput 2-d version

Definition at line 690 of file midas_vertmap.F90.

690  real, dimension(:), intent(in) :: x !< input x coordinates
691  real, dimension(:), intent(in) :: y !< input y coordinates
692  real, dimension(size(x,1),size(y,1)), intent(inout) :: x_T !< output 2-d version
693  real, dimension(size(x,1),size(y,1)), intent(inout) :: y_T !< output 2-d version
694 
695  integer :: ni,nj,i,j
696 
697  ni=size(x,1);nj=size(y,1)
698 
699  do j=1,nj
700  x_t(:,j)=x(:)
701  enddo
702 
703  do i=1,ni
704  y_t(i,:)=y(:)
705  enddo
706 
707  return
708 

◆ smooth_heights()

subroutine midas_vertmap::smooth_heights ( real, dimension(:,:), intent(inout)  zi,
integer, dimension(size(zi,1),size(zi,2)), intent(in)  fill,
integer, dimension(size(zi,1),size(zi,2)), intent(in)  bad,
real, intent(in)  sor,
integer, intent(in)  niter,
logical, intent(in)  cyclic_x,
logical, intent(in)  tripolar_n 
)
private

Solve del2 (zi) = 0 using successive iterations with a 5 point stencil. Only points fill==1 are modified. Except where bad==1, information propagates isotropically in index space. The resulting solution in each region is an approximation to del2(zi)=0 subject to boundary conditions along the valid points curve bounding this region.

Parameters
[in,out]ziinterface positions [m] or arbitrary
[in]fillpoints to be smoothed
[in]badignore these points
[in]sorsuccessive over-relaxation coefficient (typically 0.6)
[in]nitermaximum number of iterations
[in]cyclic_xinput grid cyclic condition in the zonal direction
[in]tripolar_ntripolar Arctic fold flag

Definition at line 718 of file midas_vertmap.F90.

718  real, dimension(:,:), intent(inout) :: zi !< interface positions [m] or arbitrary
719  integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill !< points to be smoothed
720  integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad !< ignore these points
721  real, intent(in) :: sor !< successive over-relaxation coefficient (typically 0.6)
722  integer, intent(in) :: niter !< maximum number of iterations
723  logical, intent(in) :: cyclic_x !< input grid cyclic condition in the zonal direction
724  logical, intent(in) :: tripolar_n !< tripolar Arctic fold flag
725 
726  integer :: i,j,k,n
727  integer :: ni,nj
728 
729  real, dimension(size(zi,1),size(zi,2)) :: res, m
730  integer, dimension(size(zi,1),size(zi,2),4) :: B
731  real, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp
732  integer, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm
733 
734  real :: Isum, bsum
735 
736  ni=size(zi,1); nj=size(zi,2)
737 
738 
739  mp=fill_boundaries(zi,cyclic_x,tripolar_n)
740 
741  b(:,:,:)=0.0
742  nm=fill_boundaries(bad,cyclic_x,tripolar_n)
743 
744  do j=1,nj
745  do i=1,ni
746  if (fill(i,j) == 1) then
747  b(i,j,1)=1-nm(i+1,j);b(i,j,2)=1-nm(i-1,j)
748  b(i,j,3)=1-nm(i,j+1);b(i,j,4)=1-nm(i,j-1)
749  endif
750  enddo
751  enddo
752 
753  do n=1,niter
754  do j=1,nj
755  do i=1,ni
756  if (fill(i,j) == 1) then
757  bsum = real(b(i,j,1)+b(i,j,2)+b(i,j,3)+b(i,j,4))
758  isum = 1.0/bsum
759  res(i,j)=isum*(b(i,j,1)*mp(i+1,j)+b(i,j,2)*mp(i-1,j)+&
760  b(i,j,3)*mp(i,j+1)+b(i,j,4)*mp(i,j-1)) - mp(i,j)
761  endif
762  enddo
763  enddo
764  res(:,:)=res(:,:)*sor
765 
766  do j=1,nj
767  do i=1,ni
768  mp(i,j)=mp(i,j)+res(i,j)
769  enddo
770  enddo
771 
772  zi(:,:)=mp(1:ni,1:nj)
773  mp = fill_boundaries(zi,cyclic_x,tripolar_n)
774  enddo
775 
776  return
777 

◆ tracer_z_init()

real function, dimension(size(tr_in,1),size(tr_in,2),nlay), public midas_vertmap::tracer_z_init ( real, dimension(:,:,:), intent(in)  tr_in,
real, dimension(size(tr_in,3)+1), intent(in)  z_edges,
real, dimension(size(tr_in,1),size(tr_in,2),nlay+1), intent(in)  e,
integer, intent(in)  nkml,
integer, intent(in)  nkbl,
real, intent(in)  land_fill,
real, dimension(size(tr_in,1),size(tr_in,2)), intent(in)  wet,
integer, intent(in)  nlay,
real, dimension(size(tr_in,1),size(tr_in,2)), intent(in), optional  nlevs,
logical, intent(in), optional  debug,
integer, intent(in), optional  i_debug,
integer, intent(in), optional  j_debug,
real, intent(in), optional  eps_z 
)

Layer model routine for remapping tracers.

Parameters
[in]tr_inThe z-space array of tracer concentrations that is read in.
[in]z_edgesThe depths of the cell edges in the input z* data [Z ~> m or m]
[in]nlayThe number of vertical layers in the target grid
[in]eThe depths of the target layer interfaces [Z ~> m or m]
[in]nkmlThe number of mixed layers
[in]nkblThe number of buffer layers
[in]land_fillfill in data over land (1)
[in]wetThe wet mask for the source data (valid points)
[in]nlevsThe number of input levels with valid data
[in]debugoptional debug flag
[in]i_debugi-index of point for debugging
[in]j_debugj-index of point for debugging
[in]eps_zA negligibly small layer thickness [Z ~> m or m].
Returns
tracers in layer space

Definition at line 153 of file midas_vertmap.F90.

153  real, dimension(:,:,:), intent(in) :: tr_in !< The z-space array of tracer concentrations that is read in.
154  real, dimension(size(tr_in,3)+1), intent(in) :: z_edges !< The depths of the cell edges in the input z* data
155  !! [Z ~> m or m]
156  integer, intent(in) :: nlay !< The number of vertical layers in the target grid
157  real, dimension(size(tr_in,1),size(tr_in,2),nlay+1), &
158  intent(in) :: e !< The depths of the target layer interfaces [Z ~> m or m]
159  integer, intent(in) :: nkml !< The number of mixed layers
160  integer, intent(in) :: nkbl !< The number of buffer layers
161  real, intent(in) :: land_fill !< fill in data over land (1)
162  real, dimension(size(tr_in,1),size(tr_in,2)), &
163  intent(in) :: wet !< The wet mask for the source data (valid points)
164  real, dimension(size(tr_in,1),size(tr_in,2)), &
165  optional, intent(in) :: nlevs !< The number of input levels with valid data
166  logical, optional, intent(in) :: debug !< optional debug flag
167  integer, optional, intent(in) :: i_debug !< i-index of point for debugging
168  integer, optional, intent(in) :: j_debug !< j-index of point for debugging
169  real, optional, intent(in) :: eps_z !< A negligibly small layer thickness [Z ~> m or m].
170  real, dimension(size(tr_in,1),size(tr_in,2),nlay) :: tr !< tracers in layer space
171 
172  ! Local variables
173  real, dimension(size(tr_in,3)) :: tr_1d !< a copy of the input tracer concentrations in a column.
174  real, dimension(nlay+1) :: e_1d ! A 1-d column of intreface heights, in the same units as e.
175  real, dimension(nlay) :: tr_ ! A 1-d column of tracer concentrations
176  integer, dimension(size(tr_in,1),size(tr_in,2)) :: nlevs_data !< number of valid levels in the input dataset
177  integer :: n,i,j,k,l,nx,ny,nz,nt,kz
178  integer :: k_top,k_bot,k_bot_prev,kk,kstart
179  real :: sl_tr ! The tracer concentration slope times the layer thickness, in tracer units.
180  real :: epsln_Z ! A negligibly thin layer thickness [Z ~> m].
181  real, dimension(size(tr_in,3)) :: wt !< The fractional weight for each layer in the range between z1 and z2
182  real, dimension(size(tr_in,3)) :: z1, z2 ! z1 and z2 are the fractional depths of the top and bottom
183  ! limits of the part of a z-cell that contributes to a layer, relative
184  ! to the cell center and normalized by the cell thickness [nondim].
185  ! Note that -1/2 <= z1 <= z2 <= 1/2.
186 
187  logical :: debug_msg, debug_, debug_pt
188 
189  nx = size(tr_in,1); ny=size(tr_in,2); nz = size(tr_in,3)
190 
191  nlevs_data = size(tr_in,3)
192  if (PRESENT(nlevs)) nlevs_data = anint(nlevs)
193  epsln_z = 1.0e-10 ; if (PRESENT(eps_z)) epsln_z = eps_z
194 
195  debug_=.false. ; if (PRESENT(debug)) debug_ = debug
196  debug_msg = debug_
197  debug_pt = debug_ ; if (PRESENT(i_debug) .and. PRESENT(j_debug)) debug_pt = debug_
198 
199  do j=1,ny
200  i_loop: do i=1,nx
201  if (nlevs_data(i,j) == 0 .or. wet(i,j) == 0.) then
202  tr(i,j,:) = land_fill
203  cycle i_loop
204  endif
205 
206  do k=1,nz
207  tr_1d(k) = tr_in(i,j,k)
208  enddo
209 
210  do k=1,nlay+1
211  e_1d(k) = e(i,j,k)
212  enddo
213  k_bot = 1 ; k_bot_prev = -1
214  do k=1,nlay
215  if (e_1d(k+1) > z_edges(1)) then
216  tr(i,j,k) = tr_1d(1)
217  elseif (e_1d(k) < z_edges(nlevs_data(i,j)+1)) then
218  if (debug_msg) then
219  print *,'*** WARNING : Found interface below valid range of z data '
220  print *,'(i,j,z_bottom,interface)= ',&
221  i,j,z_edges(nlevs_data(i,j)+1),e_1d(k)
222  print *,'z_edges= ',z_edges
223  print *,'e=',e_1d
224  print *,'*** I will extrapolate below using the bottom-most valid values'
225  debug_msg = .false.
226  endif
227  tr(i,j,k) = tr_1d(nlevs_data(i,j))
228 
229  else
230  kstart=k_bot
231  call find_overlap(z_edges, e_1d(k), e_1d(k+1), nlevs_data(i,j), &
232  kstart, k_top, k_bot, wt, z1, z2)
233 
234  if (debug_pt) then ; if ((i == i_debug) .and. (j == j_debug)) then
235  print *,'0001 k,k_top,k_bot,sum(wt),sum(z2-z1) = ',k,k_top,k_bot,sum(wt),sum(z2-z1)
236  endif ; endif
237  kz = k_top
238  sl_tr=0.0; ! cur_tr=0.0
239  if (kz /= k_bot_prev) then
240  ! Calculate the intra-cell profile.
241  if ((kz < nlevs_data(i,j)) .and. (kz > 1)) then
242  sl_tr = find_limited_slope(tr_1d, z_edges, kz)
243  endif
244  endif
245  if (kz > nlevs_data(i,j)) kz = nlevs_data(i,j)
246  ! This is the piecewise linear form.
247  tr(i,j,k) = wt(kz) * (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
248  ! For the piecewise parabolic form add the following...
249  ! + C1_3*wt(kz) * cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
250  if (debug_pt) then ; if ((i == i_debug) .and. (j == j_debug)) then
251  print *,'0002 k,k_top,k_bot,k_bot_prev,sl_tr = ',k,k_top,k_bot,k_bot_prev,sl_tr
252  endif ; endif
253 
254  do kz=k_top+1,k_bot-1
255  tr(i,j,k) = tr(i,j,k) + wt(kz)*tr_1d(kz)
256  enddo
257 
258  if (debug_pt) then ; if ((i == i_debug) .and. (j == j_debug)) then
259  print *,'0003 k,tr = ',k,tr(i,j,k)
260  endif ; endif
261 
262  if (k_bot > k_top) then
263  kz = k_bot
264  ! Calculate the intra-cell profile.
265  sl_tr = 0.0 ! ; cur_tr = 0.0
266  if ((kz < nlevs_data(i,j)) .and. (kz > 1)) then
267  sl_tr = find_limited_slope(tr_1d, z_edges, kz)
268  endif
269  ! This is the piecewise linear form.
270  tr(i,j,k) = tr(i,j,k) + wt(kz) * &
271  (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
272  ! For the piecewise parabolic form add the following...
273  ! + C1_3*cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
274 
275  if (debug_pt) then ; if ((i == i_debug) .and. (j == j_debug)) then
276  print *,'0004 k,kz,nlevs,sl_tr,tr = ',k,kz,nlevs_data(i,j),sl_tr,tr(i,j,k)
277  print *,'0005 k,kz,tr(kz),tr(kz-1),tr(kz+1) = ',k,kz,tr_1d(kz),tr_1d(kz-1),tr_1d(kz+1),z_edges(kz+2)
278  endif ; endif
279 
280  endif
281  k_bot_prev = k_bot
282 
283  endif
284  enddo ! k-loop
285 
286  do k=2,nlay ! simply fill vanished layers with adjacent value
287  if (e_1d(k)-e_1d(k+1) <= epsln_z) tr(i,j,k)=tr(i,j,k-1)
288  enddo
289 
290  enddo i_loop
291  enddo
292 

References find_limited_slope(), and find_overlap().

Here is the call graph for this function: