MOM6
Data Types | Functions/Subroutines | Variables
midas_vertmap Module Reference

Data Types

interface  fill_boundaries
 

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)
 
integer function, dimension(size(a, 1), size(x, 1)) bisect_fast (a, x, lo, hi)
 
subroutine, public determine_temperature (temp, salt, R, p_ref, niter, land_fill, h, k_start, eos)
 
subroutine find_overlap (e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z2)
 
real function find_limited_slope (val, e, k)
 
real function, dimension(size(rho, 1), size(rho, 2), size(rb, 1)), public find_interfaces (rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug)
 
subroutine, public meshgrid (x, y, x_T, y_T)
 
subroutine smooth_heights (zi, fill, bad, sor, niter, cyclic_x, tripolar_n)
 
integer function, dimension(0:size(m, 1)+1, 0:size(m, 2)+1) fill_boundaries_int (m, cyclic_x, tripolar_n)
 
real function, dimension(0:size(m, 1)+1, 0:size(m, 2)+1) fill_boundaries_real (m, cyclic_x, tripolar_n)
 

Variables

real, parameter epsln =1.e-10
 

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

Definition at line 401 of file midas_vertmap.F90.

Referenced by find_interfaces().

401 !
402 ! Return the index where to insert item x in list a, assuming a is sorted.
403 ! The return values [i] is such that all e in a[:i-1] have e <= x, and all e in
404 ! a[i:] have e > x. So if x already appears in the list, will
405 ! insert just after the rightmost x already there.
406 ! Optional args lo (default 1) and hi (default len(a)) bound the
407 ! slice of a to be searched.
408 !
409 ! (in) a - sorted list
410 ! (in) x - item to be inserted
411 ! (in) lo, hi - optional range to search
412 
413 real, dimension(:,:), intent(in) :: a
414 real, dimension(:), intent(in) :: x
415 integer, dimension(size(a,1)), intent(in), optional :: lo,hi
416 integer, dimension(size(a,1),size(x,1)) :: bi_r
417 
418 integer :: mid,num_x,num_a,i
419 integer, dimension(size(a,1)) :: lo_,hi_,lo0,hi0
420 integer :: nprofs,j
421 
422 lo_=1;hi_=size(a,2);num_x=size(x,1);bi_r=-1;nprofs=size(a,1)
423 
424 if (PRESENT(lo)) then
425  where (lo>0) lo_=lo
426 end if
427 if (PRESENT(hi)) then
428  where (hi>0) hi_=hi
429 endif
430 
431 lo0=lo_;hi0=hi_
432 
433 do j=1,nprofs
434  do i=1,num_x
435  lo_=lo0;hi_=hi0
436  do while (lo_(j) < hi_(j))
437  mid = (lo_(j)+hi_(j))/2
438  if (x(i) < a(j,mid)) then
439  hi_(j) = mid
440  else
441  lo_(j) = mid+1
442  endif
443  enddo
444  bi_r(j,i)=lo_(j)
445  enddo
446 enddo
447 
448 
449 return
450 
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), intent(in), pointer  eos 
)

Definition at line 491 of file midas_vertmap.F90.

References mom_eos::calculate_density_derivs().

Referenced by mom_state_initialization::mom_temp_salt_initialize_from_z().

491 
492 ! # This subroutine determines the potential temperature and
493 ! # salinity that is consistent with the target density
494 ! # using provided initial guess
495 ! # (inout) temp - potential temperature (degC)
496 ! # (inout) salt - salinity (PSU)
497 ! # (in) R - Desired potential density, in kg m-3.
498 ! # (in) p_ref - Reference pressure, in Pa.
499 ! # (in) niter - maximum number of iterations
500 ! # (in) h - layer thickness . Do not iterate for massless layers
501 ! # (in) k_start - starting index (i.e. below the buffer layer)
502 ! # (in) land_fill - land fill value
503 ! # (in) eos - seawater equation of state
504 
505 real, dimension(:,:,:), intent(inout) :: temp,salt
506 real, dimension(size(temp,3)), intent(in) :: r
507 real, intent(in) :: p_ref
508 integer, intent(in) :: niter
509 integer, intent(in) :: k_start
510 real, intent(in) :: land_fill
511 real, dimension(:,:,:), intent(in) :: h
512 type(eos_type), pointer, intent(in) :: eos
513 
514 real(kind=8), dimension(size(temp,1),size(temp,3)) :: t,s,dt,ds,rho,hin
515 real(kind=8), dimension(size(temp,1),size(temp,3)) :: drho_dt,drho_ds
516 real(kind=8), dimension(size(temp,1)) :: press
517 
518 integer :: nx,ny,nz,nt,i,j,k,n,itt
519 real :: dt_ds
520 logical :: adjust_salt , old_fit
521 real, parameter :: t_max = 31.0, t_min = -2.0
522 real, parameter :: s_min = 0.5, s_max=65.0
523 real, parameter :: tol=1.e-4, max_t_adj=1.0, max_s_adj = 0.5
524 
525 
526 #endif
527 
528 
529 old_fit = .true. ! reproduces siena behavior
530  ! will switch to the newer
531  ! method which simultaneously adjusts
532  ! temp and salt based on the ratio
533  ! of the thermal and haline coefficients.
534 
535 nx=size(temp,1);ny=size(temp,2); nz=size(temp,3)
536 
537 press(:) = p_ref
538 
539 do j=1,ny
540  ds(:,:) = 0. ! Needs to be zero everywhere since there is a maxval(abs(dS)) later...
541  t=temp(:,j,:)
542  s=salt(:,j,:)
543  hin=h(:,j,:)
544  dt=0.0
545  adjust_salt = .true.
546  iter_loop: do itt = 1,niter
547 #ifdef PY_SOLO
548  rho=wright_eos_2d(t,s,p_ref)
549  drho_dt=alpha_wright_eos_2d(t,s,p_ref)
550 #else
551  do k=1, nz
552  call calculate_density(t(:,k),s(:,k),press,rho(:,k),1,nx,eos)
553  call calculate_density_derivs(t(:,k),s(:,k),press,drho_dt(:,k),drho_ds(:,k),1,nx,eos)
554  enddo
555 #endif
556  do k=k_start,nz
557  do i=1,nx
558 
559 ! if (abs(rho(i,k)-R(k))>tol .and. hin(i,k)>epsln .and. abs(T(i,k)-land_fill) < epsln) then
560  if (abs(rho(i,k)-r(k))>tol) then
561  if (old_fit) then
562  dt(i,k)=(r(k)-rho(i,k))/drho_dt(i,k)
563  if (dt(i,k)>max_t_adj) dt(i,k)=max_t_adj
564  if (dt(i,k)<-1.0*max_t_adj) dt(i,k)=-1.0*max_t_adj
565  t(i,k)=max(min(t(i,k)+dt(i,k),t_max),t_min)
566  else
567  dt_ds = 10.0 - min(-drho_dt(i,k)/drho_ds(i,k),10.)
568  ds(i,k) = (r(k)-rho(i,k))/(drho_ds(i,k) - drho_dt(i,k)*dt_ds )
569  dt(i,k)= -dt_ds*ds(i,k)
570  ! if (dT(i,k)>max_t_adj) dT(i,k)=max_t_adj
571  ! if (dT(i,k)<-1.0*max_t_adj) dT(i,k)=-1.0*max_t_adj
572  t(i,k)=max(min(t(i,k)+dt(i,k),t_max),t_min)
573  s(i,k)=max(min(s(i,k)+ds(i,k),s_max),s_min)
574  endif
575  endif
576  enddo
577  enddo
578  if (maxval(abs(dt)) < tol) then
579  adjust_salt = .false.
580  exit iter_loop
581  endif
582  enddo iter_loop
583 
584  if (adjust_salt .and. old_fit) then
585  iter_loop2: do itt = 1,niter
586 #ifdef PY_SOLO
587  rho=wright_eos_2d(t,s,p_ref)
588  drho_ds=beta_wright_eos_2d(t,s,p_ref)
589 #else
590  do k=1, nz
591  call calculate_density(t(:,k),s(:,k),press,rho(:,k),1,nx,eos)
592  call calculate_density_derivs(t(:,k),s(:,k),press,drho_dt(:,k),drho_ds(:,k),1,nx,eos)
593  enddo
594 #endif
595  do k=k_start,nz
596  do i=1,nx
597 ! if (abs(rho(i,k)-R(k))>tol .and. hin(i,k)>epsln .and. abs(T(i,k)-land_fill) < epsln ) then
598  if (abs(rho(i,k)-r(k))>tol ) then
599  ds(i,k)=(r(k)-rho(i,k))/drho_ds(i,k)
600  if (ds(i,k)>max_s_adj) ds(i,k)=max_s_adj
601  if (ds(i,k)<-1.0*max_s_adj) ds(i,k)=-1.0*max_s_adj
602  s(i,k)=max(min(s(i,k)+ds(i,k),s_max),s_min)
603  endif
604  enddo
605  enddo
606  if (maxval(abs(ds)) < tol) then
607  exit iter_loop2
608  endif
609  enddo iter_loop2
610  endif
611 
612  temp(:,j,:)=t(:,:)
613  salt(:,j,:)=s(:,:)
614 enddo
615 
616 
617 return
618 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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

Definition at line 974 of file midas_vertmap.F90.

References fill_boundaries_real().

974 !
975 ! fill grid edges
976 !
977 integer, dimension(:,:), intent(in) :: m
978 logical, intent(in) :: cyclic_x, tripolar_n
979 real, dimension(size(m,1),size(m,2)) :: m_real
980 real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real
981 integer, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp
982 
983 m_real = real(m)
984 
985 mp_real = fill_boundaries_real(m_real,cyclic_x,tripolar_n)
986 
987 mp = int(mp_real)
988 
989 return
990 
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

Definition at line 994 of file midas_vertmap.F90.

Referenced by fill_boundaries_int().

994 !
995 ! fill grid edges
996 !
997 real, dimension(:,:), intent(in) :: m
998 logical, intent(in) :: cyclic_x, tripolar_n
999 real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp
1000 
1001 integer :: ni,nj,i,j
1002 
1003 ni=size(m,1); nj=size(m,2)
1004 
1005 mp(1:ni,1:nj)=m(:,:)
1006 
1007 if (cyclic_x) then
1008  mp(0,1:nj)=m(ni,1:nj)
1009  mp(ni+1,1:nj)=m(1,1:nj)
1010 else
1011  mp(0,1:nj)=m(1,1:nj)
1012  mp(ni+1,1:nj)=m(ni,1:nj)
1013 endif
1014 
1015 mp(1:ni,0)=m(1:ni,1)
1016 if (tripolar_n) then
1017  do i=1,ni
1018  mp(i,nj+1)=m(ni-i+1,nj)
1019  enddo
1020 else
1021  mp(1:ni,nj+1)=m(1:ni,nj)
1022 endif
1023 
1024 return
1025 
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 
)

Definition at line 740 of file midas_vertmap.F90.

References bisect_fast(), and epsln.

740 ! (in) rho : potential density in z-space (kg m-3)
741 ! (in) zin : levels (m)
742 ! (in) Rb : target interface densities (kg m-3)
743 ! (in) depth: ocean depth (m)
744 ! (in) nlevs: number of valid points in each column
745 ! (in) nkml : number of mixed layer pieces
746 ! (in) nkbl : number of buffer layer pieces
747 ! (in) hml : mixed layer depth
748 
749 real, dimension(:,:,:), intent(in) :: rho
750 real, dimension(size(rho,3)), intent(in) :: zin
751 real, dimension(:), intent(in) :: rb
752 real, dimension(size(rho,1),size(rho,2)), intent(in) :: depth
753 real, dimension(size(rho,1),size(rho,2)), optional, intent(in) ::nlevs
754 logical, optional, intent(in) :: debug
755 real, dimension(size(rho,1),size(rho,2),size(Rb,1)) :: zi
756 integer, intent(in), optional :: nkml, nkbl
757 real, intent(in), optional :: hml
758 
759 real, dimension(size(rho,1),size(rho,3)) :: rho_
760 real, dimension(size(rho,1)) :: depth_
761 logical :: unstable
762 integer :: dir
763 integer, dimension(size(rho,1),size(Rb,1)) :: ki_
764 real, dimension(size(rho,1),size(Rb,1)) :: zi_
765 integer, dimension(size(rho,1),size(rho,2)) :: nlevs_data
766 integer, dimension(size(rho,1)) :: lo,hi
767 real :: slope,rsm,drhodz,hml_
768 integer :: n,i,j,k,l,nx,ny,nz,nt
769 integer :: nlay,kk,nkml_,nkbl_
770 logical :: debug_ = .false.
771 
772 real, parameter :: zoff=0.999
773 
774 nlay=size(rb)-1
775 
776 zi=0.0
777 
778 
779 if (PRESENT(debug)) debug_=debug
780 
781 nx = size(rho,1); ny=size(rho,2); nz = size(rho,3)
782 nlevs_data(:,:) = size(rho,3)
783 
784 nkml_=0;nkbl_=0;hml_=0.0
785 if (PRESENT(nkml)) nkml_=max(0,nkml)
786 if (PRESENT(nkbl)) nkbl_=max(0,nkbl)
787 if (PRESENT(hml)) hml_=hml
788 
789 if (PRESENT(nlevs)) then
790  nlevs_data(:,:) = nlevs(:,:)
791 endif
792 
793 do j=1,ny
794  rho_(:,:) = rho(:,j,:)
795  i_loop: do i=1,nx
796  if (debug_) then
797  print *,'looking for interfaces, i,j,nlevs= ',i,j,nlevs_data(i,j)
798  print *,'initial density profile= ', rho_(i,:)
799  endif
800  unstable=.true.
801  dir=1
802  do while (unstable)
803  unstable=.false.
804  if (dir == 1) then
805  do k=2,nlevs_data(i,j)-1
806  if (rho_(i,k) - rho_(i,k-1) < 0.0 ) then
807  if (k.eq.2) then
808  rho_(i,k-1)=rho_(i,k)-epsln
809  else
810  drhodz = (rho_(i,k+1)-rho_(i,k-1))/(zin(k+1)-zin(k-1))
811  if (drhodz < 0.0) then
812  unstable=.true.
813  endif
814  rho_(i,k) = rho_(i,k-1)+drhodz*zoff*(zin(k)-zin(k-1))
815  endif
816  endif
817  enddo
818  dir=-1*dir
819  else
820  do k=nlevs_data(i,j)-1,2,-1
821  if (rho_(i,k+1) - rho_(i,k) < 0.0) then
822  if (k .eq. nlevs_data(i,j)-1) then
823  rho_(i,k+1)=rho_(i,k-1)+epsln
824  else
825  drhodz = (rho_(i,k+1)-rho_(i,k-1))/(zin(k+1)-zin(k-1))
826  if (drhodz < 0.0) then
827  unstable=.true.
828  endif
829  rho_(i,k) = rho_(i,k+1)-drhodz*(zin(k+1)-zin(k))
830  endif
831  endif
832  enddo
833  dir=-1*dir
834  endif
835  enddo
836  if (debug_) then
837  print *,'final density profile= ', rho_(i,:)
838  endif
839  enddo i_loop
840 
841  ki_(:,:) = 0
842  zi_(:,:) = 0.0
843  depth_(:)=-1.0*depth(:,j)
844  lo(:)=1
845  hi(:)=nlevs_data(:,j)
846  ki_ = bisect_fast(rho_,rb,lo,hi)
847  ki_(:,:) = max(1,ki_(:,:)-1)
848  do i=1,nx
849  do l=2,nlay
850  slope = (zin(ki_(i,l)+1) - zin(ki_(i,l)))/max(rho_(i,ki_(i,l)+1) - rho_(i,ki_(i,l)),epsln)
851  zi_(i,l) = -1.0*(zin(ki_(i,l)) + slope*(rb(l)-rho_(i,ki_(i,l))))
852  zi_(i,l) = max(zi_(i,l),depth_(i))
853  zi_(i,l) = min(zi_(i,l),-1.0*hml_)
854  enddo
855  zi_(i,nlay+1)=depth_(i)
856  do l=2,nkml_+1
857  zi_(i,l)=max(((1.0-real(l))/real(nkml_))*hml_,depth_(i))
858  enddo
859  do l=nlay,nkml_+2,-1
860  if (zi_(i,l) < zi_(i,l+1)+epsln) then
861  zi_(i,l)=zi_(i,l+1)+epsln
862  endif
863  if (zi_(i,l)>-1.0*hml_) then
864  zi_(i,l)=max(-1.0*hml_,depth_(i))
865  endif
866  enddo
867  enddo
868  zi(:,j,:)=zi_(:,:)
869 enddo
870 
871 return
872 
873 
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

Definition at line 694 of file midas_vertmap.F90.

Referenced by tracer_z_init().

694 
695 ! This subroutine determines a limited slope for val to be advected with
696 ! a piecewise limited scheme.
697 
698 ! Arguments: val - An column the values that are being interpolated.
699 ! (in) e - A column's interface heights, in m.
700 ! (in) slope - The normalized slope in the intracell distribution of val.
701 ! (in) k - The layer whose slope is being determined.
702 
703 
704 real, dimension(:), intent(in) :: val
705 real, dimension(:), intent(in) :: e
706 integer, intent(in) :: k
707 real :: slope,amx,bmx,amn,bmn,cmn,dmn
708 
709 real :: d1, d2
710 
711 if ((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) then
712  slope = 0.0 ! ; curvature = 0.0
713 else
714  d1 = 0.5*(e(k-1)-e(k+1)) ; d2 = 0.5*(e(k)-e(k+2))
715  slope = ((d1**2)*(val(k+1) - val(k)) + (d2**2)*(val(k) - val(k-1))) * &
716  (e(k) - e(k+1)) / (d1*d2*(d1+d2))
717 ! slope = 0.5*(val(k+1) - val(k-1))
718 ! This is S.J. Lin's form of the PLM limiter.
719  amx=max(val(k-1),val(k))
720  bmx = max(amx,val(k+1))
721  amn = min(abs(slope),2.0*(bmx-val(k)))
722  bmn = min(val(k-1),val(k))
723  cmn = 2.0*(val(k)-min(bmn,val(k+1)))
724  dmn = min(amn,cmn)
725  slope = sign(1.0,slope) * dmn
726 
727 ! min(abs(slope), &
728 ! 2.0*(max(val(k-1),val(k),val(k+1)) - val(k)), &
729 ! 2.0*(val(k) - min(val(k-1),val(k),val(k+1))))
730 ! curvature = 0.0
731 endif
732 
733 return
734 
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

Definition at line 623 of file midas_vertmap.F90.

Referenced by tracer_z_init().

623 
624 ! This subroutine determines the layers bounded by interfaces e that overlap
625 ! with the depth range between Z_top and Z_bot, and also the fractional weights
626 ! of each layer. It also calculates the normalized relative depths of the range
627 ! of each layer that overlaps that depth range.
628 ! Note that by convention, e decreases with increasing k and Z_top > Z_bot.
629 !
630 ! Arguments: e - A column's interface heights, in m.
631 ! (in) Z_top - The top of the range being mapped to, in m.
632 ! (in) Z_bot - The bottom of the range being mapped to, in m.
633 ! (in) k_max - The number of valid layers.
634 ! (in) k_start - The layer at which to start searching.
635 ! (out) k_top, k_bot - The indices of the top and bottom layers that
636 ! overlap with the depth range.
637 ! (out) wt - The relative weights of each layer from k_top to k_bot.
638 ! (out) z1, z2 - z1 and z2 are the depths of the top and bottom limits of
639 ! the part of a layer that contributes to a depth level,
640 ! relative to the cell center and normalized by the cell
641 ! thickness, nondim. Note that -1/2 <= z1 < z2 <= 1/2.
642 
643 real, dimension(:), intent(in) :: e
644 real, intent(in) :: z_top, z_bot
645 integer, intent(in) :: k_max, k_start
646 integer, intent(out) :: k_top, k_bot
647 real, dimension(:), intent(out) :: wt, z1, z2
648 
649 real :: ih, e_c, tot_wt, i_totwt
650 integer :: k
651 
652 wt(:)=0.0; z1(:)=0.0; z2(:)=0.0
653 k_top = k_start; k_bot= k_start; wt(1) = 1.0; z1(1)=-0.5; z2(1) = 0.5
654 
655 do k=k_start,k_max ;if (e(k+1)<z_top) exit ; enddo
656 k_top = k
657 if (k>k_max) return
658 
659 ! Determine the fractional weights of each layer.
660 ! Note that by convention, e and Z_int decrease with increasing k.
661 if (e(k+1)<=z_bot) then
662  wt(k) = 1.0 ; k_bot = k
663  ih = 1.0 / (e(k)-e(k+1))
664  e_c = 0.5*(e(k)+e(k+1))
665  z1(k) = (e_c - min(e(k),z_top)) * ih
666  z2(k) = (e_c - z_bot) * ih
667 else
668  wt(k) = min(e(k),z_top) - e(k+1) ; tot_wt = wt(k) ! These are always > 0.
669  z1(k) = (0.5*(e(k)+e(k+1)) - min(e(k),z_top)) / (e(k)-e(k+1))
670  z2(k) = 0.5
671  k_bot = k_max
672  do k=k_top+1,k_max
673  if (e(k+1)<=z_bot) then
674  k_bot = k
675  wt(k) = e(k) - z_bot ; z1(k) = -0.5
676  z2(k) = (0.5*(e(k)+e(k+1)) - z_bot) / (e(k)-e(k+1))
677  else
678  wt(k) = e(k) - e(k+1) ; z1(k) = -0.5 ; z2(k) = 0.5
679  endif
680  tot_wt = tot_wt + wt(k) ! wt(k) is always > 0.
681  if (k>=k_bot) exit
682  enddo
683 
684  i_totwt = 1.0 / tot_wt
685  do k=k_top,k_bot ; wt(k) = i_totwt*wt(k) ; enddo
686 endif
687 
688 return
689 
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 
)

Definition at line 877 of file midas_vertmap.F90.

877 
878 ! create a 2d-mesh of grid coordinates
879 ! from 1-d arrays.
880 
881 real, dimension(:), intent(in) :: x,y
882 real, dimension(size(x,1),size(y,1)), intent(inout) :: x_t,y_t
883 
884 integer :: ni,nj,i,j
885 
886 ni=size(x,1);nj=size(y,1)
887 
888 do j=1,nj
889  x_t(:,j)=x(:)
890 enddo
891 
892 do i=1,ni
893  y_t(i,:)=y(:)
894 enddo
895 
896 return
897 

◆ 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

Definition at line 901 of file midas_vertmap.F90.

901 !
902 ! Solve del2 (zi) = 0 using successive iterations
903 ! with a 5 point stencil. Only points fill==1 are
904 ! modified. Except where bad==1, information propagates
905 ! isotropically in index space. The resulting solution
906 ! in each region is an approximation to del2(zi)=0 subject to
907 ! boundary conditions along the valid points curve bounding this region.
908 
909 
910 real, dimension(:,:), intent(inout) :: zi
911 integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill
912 integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad
913 real, intent(in) :: sor
914 integer, intent(in) :: niter
915 logical, intent(in) :: cyclic_x, tripolar_n
916 
917 integer :: i,j,k,n
918 integer :: ni,nj
919 
920 real, dimension(size(zi,1),size(zi,2)) :: res, m
921 integer, dimension(size(zi,1),size(zi,2),4) :: b
922 real, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp
923 integer, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm
924 
925 real :: isum, bsum
926 
927 ni=size(zi,1); nj=size(zi,2)
928 
929 
930 mp=fill_boundaries(zi,cyclic_x,tripolar_n)
931 
932 b(:,:,:)=0.0
933 nm=fill_boundaries(bad,cyclic_x,tripolar_n)
934 
935 do j=1,nj
936  do i=1,ni
937  if (fill(i,j) .eq. 1) then
938  b(i,j,1)=1-nm(i+1,j);b(i,j,2)=1-nm(i-1,j)
939  b(i,j,3)=1-nm(i,j+1);b(i,j,4)=1-nm(i,j-1)
940  endif
941  enddo
942 enddo
943 
944 do n=1,niter
945  do j=1,nj
946  do i=1,ni
947  if (fill(i,j) .eq. 1) then
948  bsum = real(b(i,j,1)+b(i,j,2)+b(i,j,3)+b(i,j,4))
949  isum = 1.0/bsum
950  res(i,j)=isum*(b(i,j,1)*mp(i+1,j)+b(i,j,2)*mp(i-1,j)+&
951  b(i,j,3)*mp(i,j+1)+b(i,j,4)*mp(i,j-1)) - mp(i,j)
952  endif
953  enddo
954  enddo
955  res(:,:)=res(:,:)*sor
956 
957  do j=1,nj
958  do i=1,ni
959  mp(i,j)=mp(i,j)+res(i,j)
960  enddo
961  enddo
962 
963  zi(:,:)=mp(1:ni,1:nj)
964  mp = fill_boundaries(zi,cyclic_x,tripolar_n)
965 end do
966 
967 
968 
969 return
970 

◆ 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 
)

Definition at line 214 of file midas_vertmap.F90.

References epsln, find_limited_slope(), and find_overlap().

214 !
215 ! Adopted from R. Hallberg
216 ! Arguments:
217 ! (in) tr_in - The z-space array of tracer concentrations that is read in.
218 ! (in) z_edges - The depths of the cell edges in the input z* data (m)
219 ! (in) e - The depths of the layer interfaces (m)
220 ! (in) nkml - number of mixed layer pieces
221 ! (in) nkbl - number of buffer layer pieces
222 ! (in) land_fill - fill in data over land
223 ! (in) wet - wet mask (1=ocean)
224 ! (in) nlay - number of layers
225 ! (in) nlevs - number of levels
226 
227 ! (out) tr - tracers on layers
228 
229 ! tr_1d ! A copy of the input tracer concentrations in a column.
230 ! wt ! The fractional weight for each layer in the range between
231  ! k_top and k_bot, nondim.
232 ! z1 ! z1 and z2 are the depths of the top and bottom limits of the part
233 ! z2 ! of a z-cell that contributes to a layer, relative to the cell
234 ! center and normalized by the cell thickness, nondim.
235 ! Note that -1/2 <= z1 <= z2 <= 1/2.
236 !
237 real, dimension(:,:,:), intent(in) :: tr_in
238 real, dimension(size(tr_in,3)+1), intent(in) :: z_edges
239 integer, intent(in) :: nlay
240 real, dimension(size(tr_in,1),size(tr_in,2),nlay+1), intent(in) :: e
241 integer, intent(in) :: nkml,nkbl
242 real, intent(in) :: land_fill
243 real, dimension(size(tr_in,1),size(tr_in,2)), intent(in) :: wet
244 real, dimension(size(tr_in,1),size(tr_in,2)), optional, intent(in) ::nlevs
245 logical, intent(in), optional :: debug
246 integer, intent(in), optional :: i_debug, j_debug
247 
248 real, dimension(size(tr_in,1),size(tr_in,2),nlay) :: tr
249 real, dimension(size(tr_in,3)) :: tr_1d
250 real, dimension(nlay+1) :: e_1d
251 real, dimension(nlay) :: tr_
252 integer, dimension(size(tr_in,1),size(tr_in,2)) :: nlevs_data
253 
254 integer :: n,i,j,k,l,nx,ny,nz,nt,kz
255 integer :: k_top,k_bot,k_bot_prev,kk,kstart
256 real :: sl_tr
257 real, dimension(size(tr_in,3)) :: wt,z1,z2
258 logical :: debug_msg, debug_
259 
260 nx = size(tr_in,1); ny=size(tr_in,2); nz = size(tr_in,3)
261 
262 nlevs_data = size(tr_in,3)
263 if (PRESENT(nlevs)) then
264  nlevs_data = anint(nlevs)
265 endif
266 
267 debug_=.false.
268 if (PRESENT(debug)) then
269  debug_=debug
270 endif
271 
272 debug_msg = .false.
273 if (debug_) then
274  debug_msg=.true.
275 endif
276 
277 
278 do j=1,ny
279  i_loop: do i=1,nx
280  if (nlevs_data(i,j) .eq. 0 .or. wet(i,j) .eq. 0.) then
281  tr(i,j,:) = land_fill
282  cycle i_loop
283  endif
284 
285  do k=1,nz
286  tr_1d(k) = tr_in(i,j,k)
287  enddo
288 
289  do k=1,nlay+1
290  e_1d(k) = e(i,j,k)
291  enddo
292  k_bot = 1 ; k_bot_prev = -1
293  do k=1,nlay
294  if (e_1d(k+1) > z_edges(1)) then
295  tr(i,j,k) = tr_1d(1)
296  elseif (e_1d(k) < z_edges(nlevs_data(i,j)+1)) then
297  if (debug_msg) then
298  print *,'*** WARNING : Found interface below valid range of z data '
299  print *,'(i,j,z_bottom,interface)= ',&
300  i,j,z_edges(nlevs_data(i,j)+1),e_1d(k)
301  print *,'z_edges= ',z_edges
302  print *,'e=',e_1d
303  print *,'*** I will extrapolate below using the bottom-most valid values'
304  debug_msg = .false.
305  endif
306  tr(i,j,k) = tr_1d(nlevs_data(i,j))
307 
308  else
309  kstart=k_bot
310  call find_overlap(z_edges, e_1d(k), e_1d(k+1), nlevs_data(i,j), &
311  kstart, k_top, k_bot, wt, z1, z2)
312 
313  if (debug_) then
314  if (PRESENT(i_debug)) then
315  if (i.eq.i_debug.and.j.eq.j_debug) then
316  print *,'0001 k,k_top,k_bot,sum(wt),sum(z2-z1) = ',k,k_top,k_bot,sum(wt),sum(z2-z1)
317  endif
318  endif
319  endif
320  kz = k_top
321  sl_tr=0.0; ! cur_tr=0.0
322  if (kz /= k_bot_prev) then
323 ! Calculate the intra-cell profile.
324  if ((kz < nlevs_data(i,j)) .and. (kz > 1)) then
325  sl_tr = find_limited_slope(tr_1d, z_edges, kz)
326  endif
327  endif
328  if (kz > nlevs_data(i,j)) kz = nlevs_data(i,j)
329 ! This is the piecewise linear form.
330  tr(i,j,k) = wt(kz) * &
331  (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
332 ! For the piecewise parabolic form add the following...
333 ! + C1_3*cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
334 ! if (debug_) then
335 ! print *,'k,k_top,k_bot= ',k,k_top,k_bot
336 ! endif
337  if (debug_) then
338  if (PRESENT(i_debug)) then
339  if (i.eq.i_debug.and.j.eq.j_debug) then
340  print *,'0002 k,k_top,k_bot,k_bot_prev,sl_tr = ',k,k_top,k_bot,k_bot_prev,sl_tr
341  endif
342  endif
343  endif
344 
345  do kz=k_top+1,k_bot-1
346  tr(i,j,k) = tr(i,j,k) + wt(kz)*tr_1d(kz)
347  enddo
348 
349  if (debug_) then
350  if (PRESENT(i_debug)) then
351  if (i.eq.i_debug.and.j.eq.j_debug) then
352  print *,'0003 k,tr = ',k,tr(i,j,k)
353  endif
354  endif
355  endif
356 
357  if (k_bot > k_top) then
358  kz = k_bot
359 ! Calculate the intra-cell profile.
360  sl_tr = 0.0 ! ; cur_tr = 0.0
361  if ((kz < nlevs_data(i,j)) .and. (kz > 1)) then
362  sl_tr = find_limited_slope(tr_1d, z_edges, kz)
363 ! if (debug_) then
364 ! print *,'002 sl_tr,k,kz,nlevs= ',sl_tr,k,kz,nlevs_data(i,j),nlevs(i,j)
365 ! endif
366  endif
367 ! This is the piecewise linear form.
368  tr(i,j,k) = tr(i,j,k) + wt(kz) * &
369  (tr_1d(kz) + 0.5*sl_tr*(z2(kz) + z1(kz)))
370 ! For the piecewise parabolic form add the following...
371 ! + C1_3*cur_tr*(z2(kz)**2 + z2(kz)*z1(kz) + z1(kz)**2))
372 
373  if (debug_) then
374  if (PRESENT(i_debug)) then
375  if (i.eq.i_debug.and.j.eq.j_debug) then
376  print *,'0004 k,kz,nlevs,sl_tr,tr = ',k,kz,nlevs_data(i,j),sl_tr,tr(i,j,k)
377  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)
378  endif
379  endif
380  endif
381 
382  endif
383  k_bot_prev = k_bot
384 
385  endif
386  enddo ! k-loop
387 
388  do k=2,nlay ! simply fill vanished layers with adjacent value
389  if (e_1d(k)-e_1d(k+1) .le. epsln) tr(i,j,k)=tr(i,j,k-1)
390  enddo
391 
392  enddo i_loop
393 enddo
394 
395 return
396 
Here is the call graph for this function:

Variable Documentation

◆ epsln

real, parameter midas_vertmap::epsln =1.e-10
private

Definition at line 63 of file midas_vertmap.F90.

Referenced by find_interfaces(), and tracer_z_init().

63  real, parameter :: epsln=1.e-10