module hmix_GM !BOP ! !MODULE: hmix_gm ! !DESCRIPTION: ! This module contains routines for computing horizontal mixing ! using the Gent-McWilliams eddy transport parameterization ! and isopycnal diffusion. ! !REVISION HISTORY: ! 2018-09-04 ! !USES: use const use physics use hmix_GM_submeso_share implicit none private ! save ! !PUBLIC MEMBER FUNCTIONS: public :: init_gm, & hdifft_gm !EOC !BOC !----------------------------------------------------------------------- ! ! variables to save from one call to next ! !----------------------------------------------------------------------- integer, parameter :: & ktp = 1, kbt = 2 ! refer to the top and bottom halves of a ! grid cell, respectively real*8, dimension(i0,j0) :: & HYXW, HXYS, &! west and south-shifted values of above RB, &! Rossby radius RBR, &! inverse of Rossby radius BL_DEPTH, &! boundary layer depth UIT, VIT, &! work arrays for isopycnal mixing velocities WTOP_ISOP, WBOT_ISOP ! vertical component of isopycnal velocities real*8, dimension(i0,j0,2,2,k1) :: & SF_SLX, SF_SLY ! components of the merged streamfunction real*8, dimension(i0,j0,2,k1) :: & SLA_SAVE ! isopycnal slopes real*8, dimension(i0,j0,tnt) :: & FZTOP ! vertical flux logical :: & cancellation_occurs ! specified choices for the isopycnal and ! thickness diffusion coefficients ! result in cancellation of some tensor ! elements real*8, dimension(0:k1) :: dzw, dzwr !----------------------------------------------------------------------- ! ! KAPPA_LATERAL and KAPPA_VERTICAL are 2D and 3D arrays, ! respectively, ! containing the spatial variations of the isopycnal (KAPPA_ISOP) ! and thickness (KAPPA_THIC) diffusion coefficients. Except in ! kappa_type_eg, ! KAPPA_LATERAL has the actual diffusion coefficient values in ! cm^2/s, ! whereas KAPPA_VERTICAL is unitless. So, the total coefficients are ! constructed using ! ! KAPPA_ISOP or KAPPA_THIC (:,:,:,k,bid) ~ KAPPA_LATERAL(:,:,bid) ! * ! KAPPA_VERTICAL(:,:,k,bid) ! ! When kappa_type_eg, KAPPA_VERTICAL contains the isopycnal ! diffusivity ! coefficients in cm^2/s and KAPPA_LATERAL is not used! ! !----------------------------------------------------------------------- ! real*8, dimension(i0,j0,2,k1) :: & ! KAPPA_ISOP, & ! 3D isopycnal diffusion coefficient ! ! for top and bottom half of a grid cell ! KAPPA_THIC, & ! 3D thickness diffusion coefficient ! ! for top and bottom half of a grid cell ! HOR_DIFF ! 3D horizontal diffusion coefficient ! ! for top and bottom half of a grid cell real*8, dimension(i0,j0) :: & KAPPA_LATERAL ! horizontal variation of KAPPA in cm^2/s real*8, dimension(i0,j0,k1) :: & KAPPA_VERTICAL ! vertical variation of KAPPA (unitless), ! e.g. normalized buoyancy frequency ! dependent ! profiles at the tracer grid points ! ( = N^2 / N_ref^2 ) OR a time-independent ! user-specified function ! real*8, dimension(i0,j0,k1) :: VDC_GM real*8, dimension(i0,j0,k1) :: & BUOY_FREQ_SQ ! N^2 defined at level interfaces ! true for diagnostic bolus velocity computation logical, parameter :: diag_gm_bolus = .true. !----------------------------------------------------------------------- ! ! if use_const_ah_bkg_srfbl = .true., then the specified constant ! value of ah_bkg_srfbl is used as the "maximum" background ! horizontal ! diffusivity within the surface boundary layer. Otherwise, ! KAPPA_ISOP is utilized as this "maximum". ! !----------------------------------------------------------------------- logical, parameter :: use_const_ah_bkg_srfbl = .true. ! control for transition layer parameterization logical, parameter :: transition_layer_on = .true. ! isopycnal diffusivity real*8, parameter :: ah = 3.0e7 ! thickness (GM bolus) diffusivity real*8, parameter :: ah_bolus = 3.0e7 ! backgroud horizontal diffusivity at k= KMT real*8, parameter :: ah_bkg_bottom = 0.d0 ! backgroud horizontal diffusivity within the surface boundary layer real*8, parameter :: ah_bkg_srfbl = 3.0e7 ! max. slope allowed for redi diffusion ! real*8, parameter :: slm_r = 0.3d0 real*8, parameter :: slm_r = 0.3d0 ! max. slope allowed for bolus transport ! real*8, parameter :: slm_b = 0.3d0 real*8, parameter :: slm_b = 0.3d0 !----------------------------------------------------------------------- ! ! transition layer type variables ! !!----------------------------------------------------------------------- type tlt_info real*8, dimension(i0,j0) :: & DIABATIC_DEPTH, & ! depth of the diabatic region at the ! surface, i.e. mean mixed or boundary ! layer depth THICKNESS, & ! transition layer thickness INTERIOR_DEPTH ! depth at which the interior, adiabatic ! region starts, i.e. ! = TLT%DIABATIC_DEPTH + TLT%THICKNESS integer, dimension(i0,j0) :: & K_LEVEL, & ! k level at or below which the interior, ! adiabatic region starts ZTW ! designates if the interior region ! starts below depth zt or zw. ! ( = 1 for zt, = 2 for zw ) end type tlt_info type (tlt_info) :: & TLT ! transition layer thickness related fields !*********************************************************************** contains !*********************************************************************** !BOP ! !IROUTINE: init_gm ! !INTERFACE: subroutine init_gm ! !DESCRIPTION: ! Initializes various choices and allocates necessary space. ! Also computes some time-independent arrays. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! input namelist for setting various GM options ! !----------------------------------------------------------------------- integer :: & k, &! vertical level index i, j ! lateral indices ! real*8, parameter :: & ! ah = 0.8e7, & ! ah_bolus = 0.8e7, &!thickness diffusion coefficient(cm^2/s) ! ah_bkg_bottom = 0.d0, &!background horizontal diffusion at k=KMT ! ah_bkg_srfbl = 0.8e7, &!background horizontal diffusion within the !surface boundary layer (cm^2/s) ! slm_r = 0.3d0, &!max. slope allowed for isopycnal(redi) diffusion (= 0.3) ! slm_b = 0.3d0 !max. slope allowed for thickness(bolus) diffusion (= 0.3) !for diagnostic bolus velocity computation ! logical, parameter :: diag_gm_bolus = .true. ! logical, parameter :: use_const_ah_bkg_srfbl = .true. ! logical, parameter :: transition_layer_on = .true. real*8, dimension(i0,j0) :: FCORT0 ! coriolis parameter !----------------------------------------------------------------------- ! ! initialize bolus velocity arrays if requested ! !----------------------------------------------------------------------- HYXW = 0.d0 HXYS = 0.d0 RBR = 0.d0 BL_DEPTH = 0.d0 SF_SLX = 0.d0 SF_SLY = 0.d0 FZTOP = 0.d0 VDC_GM = 0.d0 SLA_SAVE = 0.d0 RB = 0.d0 WTOP_ISOP = 0.d0 WBOT_ISOP = 0.d0 UIT = 0.d0 VIT = 0.d0 !----------------------------------------------------------------------- ! ! initialize various time-independent arrays ! !----------------------------------------------------------------------- KAPPA_LATERAL(:,:) = ah KAPPA_VERTICAL(:,:,:) = 1.0d0 KAPPA_ISOP(:,:,:,:) = ah KAPPA_THIC(:,:,:,:) = ah_bolus HOR_DIFF(:,:,:,:) = ah HYXW(:,:) = eoshift(HYX(:,:), dim=1, shift=-1) HXYS(:,:) = eoshift(HXY(:,:), dim=2, shift=-1) dzw(0) = 0.5d0*(z(3) - z(1)) dzw(k1) = 0.5d0*(z(2*k1+1) - z(2*k1-1)) dzwr(0) = 1.0d0/dzw(0) dzwr(k1) = 1.0d0/dzw(k1) do k = 1,k1-1 !dzw(k) = 1.0d0/odzw(k+1) dzw(k) = z(2*k+2) - z(2*k) !dzwr(k) = odzw(k+1) dzwr(k) = 1.0d0/dzw(k) enddo !----------------------------------------------------------------------- ! compute the Rossby radius which will be used to ! contol KAPPA [Large et al (1997), JPO, 27, ! pp 2418-2447]. Rossby radius = c/(2*omg*sin(latitude)) ! where c=200cm/s is the first baroclinic wave speed. ! 15km < Rossby radius < 100km !----------------------------------------------------------------------- !*** Inverse of Rossby radius do j=1,j0 do i=1,i0 FCORT0(i,j) = fcort(j) enddo enddo RBR(:,:) = abs(FCORT0(:,:))/200.0 ! |f|/Cg, Cg = 2 m/s = 200 cm/s RBR(:,:) = min(RBR(:,:), 1.0d0/1.5e+6)! Cg/|f| .ge. 15 km = 1.5e+6 cm RBR(:,:) = max(RBR(:,:),1.e-7) ! Cg/|f| .le. 100 km = 1.e+7 cm if ( transition_layer_on ) then RB(:,:) = 1.0d0 / RBR(:,:) endif BUOY_FREQ_SQ = 0.d0 !*** for transition layer cases, the following will always be true!! if ( transition_layer_on ) cancellation_occurs = .false. !----------------------------------------------------------------------- !EOC end subroutine init_gm !*********************************************************************** !BOP ! !IROUTINE: hdifft_gm ! !INTERFACE: subroutine hdifft_gm (k, TMIX, vdc, HOR_DIFF_X, HOR_DIFF_Y, & KH_GM_X, KH_GM_Y, GTK0) ! !DESCRIPTION: ! Gent-McWilliams eddy transport parameterization ! and isopycnal diffusion. ! ! This routine must be called successively with k = 1,2,3,... ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: integer, intent(in) :: k ! depth level index real*8, dimension(i0,j0,k1,tnt), intent(in) :: & TMIX ! tracers at all vertical levels ! at mixing time level real*8, dimension(i0,j0,0:k0,2),intent(inout) :: & vdc ! diffusivity for tracer diffusion ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0), intent(out) :: & HOR_DIFF_X, & HOR_DIFF_Y, & KH_GM_X, &! isopycanl + background horrizontal diffusion KH_GM_y ! coefficients used for calculating horizontal flux ! in x- or y-direction thru vertical faces real*8, dimension(i0,j0,tnt), intent(out) :: & ! GTK &! diffusion+bolus advection for nth tracer at level k GTK0 ! diffusion+bolus advection added in MUSOC_v1.1 !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer, parameter :: & ieast = 1, iwest = 2, & jnorth = 1, jsouth = 2 integer :: & i,j,n,kk, &! dummy loop counters kid, &! array indices kk_sub, kp1 real*8 :: & fz, dz_bottom, factor real*8, dimension(i0,j0) :: & CX, CY, & RZ, &! Dz(rho) SLA, &! absolute value of slope WORK1, WORK2, &! local work space WORK3, WORK4, &! local work space KMASK, &! ocean mask TAPER1, TAPER2, TAPER3, &! tapering factors UIB, VIB, &! work arrays for isopycnal mixing velocities U_ISOP, V_ISOP ! horizontal components of isopycnal velocities real*8, dimension(i0,j0,tnt) :: & FX, FY, FX0, FY0 ! fluxes across east, north faces real*8, dimension(2) :: & reference_depth !----------------------------------------------------------------------- ! ! initialize various quantities ! !----------------------------------------------------------------------- !if (myid==0) write(*,*) 'check=', transition_layer_on, diag_gm_bolus !if (myid==0) write(*,*) 'check=', use_const_ah_bkg_srfbl, cancellation_occurs U_ISOP = 0.d0 V_ISOP = 0.d0 WORK1 = 0.d0 WORK2 = 0.d0 WORK3 = 0.d0 WORK4 = 0.d0 if ( k == 1 ) then if ( diag_gm_bolus ) then UIB = 0.d0 VIB = 0.d0 UIT(:,:) = 0.d0 VIT(:,:) = 0.d0 WBOT_ISOP(:,:) = 0.d0 endif HOR_DIFF(:,:,ktp,k) = ah_bkg_srfbl BL_DEPTH(:,:) = kpp_hblt(:,:) endif if (diag_gm_bolus) WTOP_ISOP(:,:) = WBOT_ISOP(:,:) ! CX = merge(HYX(:,:)*0.25d0, 0.d0, (k <= KMT(:,:)) & ! .and. (k <= KMTE(:,:))) ! CY = merge(HXY(:,:)*0.25d0, 0.d0, (k <= KMT(:,:)) & ! .and. (k <= KMTN(:,:))) if (k==1) then CX = HYX*0.25d0*iu0(1:i0,:) CY = HXY*0.25d0*iv0(:,1:j0) else CX = HYX*0.25d0*iu(1:i0,:,k) CY = HXY*0.25d0*iv(:,1:j0,k) endif if ( k == 1 ) then if ( transition_layer_on ) then call smooth_hblt ( .false., .true., & smooth_out=TLT%DIABATIC_DEPTH(:,:) ) do kk=1,k1 do kk_sub = ktp,kbt kid = kk + kk_sub - 2 SLA_SAVE(:,:,kk_sub,kk) = dzw(kid)*sqrt(.5d0*( & (SLX(:,:,1,kk_sub,kk)**2 & + SLX(:,:,2,kk_sub,kk)**2)/DXT(:,:)**2 & + (SLY(:,:,1,kk_sub,kk)**2 & + SLY(:,:,2,kk_sub,kk)**2)/DYT(:,:)**2)) & + eps enddo enddo call transition_layer endif !----------------------------------------------------------------------- ! ! compute isopycnal and thickness diffusion coefficients if ! they depend on the model fields ! !----------------------------------------------------------------------- call buoyancy_frequency_dependent_profile (TMIX) !----------------------------------------------------------------------- ! ! reinitialize the diffusivity coefficients ! !----------------------------------------------------------------------- do kk_sub=ktp,kbt do kk=1,k1 KAPPA_ISOP(:,:,kk_sub,kk) = KAPPA_LATERAL(:,:) & * KAPPA_VERTICAL(:,:,kk) enddo enddo do kk_sub=ktp,kbt do kk=1,k1 KAPPA_THIC(:,:,kk_sub,kk) = ah_bolus & * KAPPA_VERTICAL(:,:,kk) enddo enddo !----------------------------------------------------------------------- ! ! control slope of isopycnal surfaces or KAPPA ! !----------------------------------------------------------------------- do kk=1,k1 kp1 = min(kk+1,k1) reference_depth(ktp) = z(2*kp1) -z(1) reference_depth(kbt) = z(2*kp1+1) - z(1) if ( kk == k1 ) reference_depth(ktp) = z(2*kp1+1) - z(1) do kk_sub = ktp,kbt kid = kk + kk_sub - 2 !----------------------------------------------------------------------- ! ! control KAPPA to reduce the isopycnal mixing near the ! ocean surface Large et al (1997), JPO, 27, pp 2418-2447. ! WORK1 = ratio between the depth of water parcel and ! the vertical displacement of isopycnal surfaces ! where the vertical displacement = ! Rossby radius * slope of isopycnal surfaces ! !----------------------------------------------------------------------- if ( transition_layer_on ) then SLA = SLA_SAVE(:,:,kk_sub,kk) endif TAPER1 = 1.0d0 !----------------------------------------------------------------------- ! ! control KAPPA for numerical stability ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! methods to control slope ! case (slope_control_notanh) ! diff_tapering = .false. !----------------------------------------------------------------------- TAPER2 = 1.0d0 TAPER3 = 1.0d0 ! similar to DM95 except replacing tanh by ! function = x*(1.-0.25*abs(x)) for |x|<2 ! = sign(x) for |x|>2 ! (faster than DM95) do j=1,j0 do i=1,i0 if (SLA(i,j) > 0.2d0*slm_r .and. & SLA(i,j) < 0.6d0*slm_r) then TAPER2(i,j) = & 0.5d0*(1.0d0-(2.5d0*SLA(i,j)/slm_r-1.0d0)* & (4.0d0-abs(10.0d0*SLA(i,j)/slm_r-4.0d0))) else if (SLA(i,j) >= 0.6d0*slm_r) then TAPER2(i,j) = 0.d0 endif enddo enddo TAPER3 = TAPER2 if ( transition_layer_on ) then TAPER2 = merge(1.0d0, TAPER2, reference_depth(kk_sub) & <= TLT%DIABATIC_DEPTH(:,:)) TAPER3 = merge(1.0d0, TAPER3, reference_depth(kk_sub) & <= TLT%DIABATIC_DEPTH(:,:)) endif if ( transition_layer_on .and. use_const_ah_bkg_srfbl ) then HOR_DIFF(:,:,kk_sub,kk) = ah_bkg_srfbl endif KAPPA_ISOP(:,:,kk_sub,kk) = & TAPER1 * TAPER2 * KAPPA_ISOP(:,:,kk_sub,kk) KAPPA_THIC(:,:,kk_sub,kk) = & TAPER1 * TAPER3 * KAPPA_THIC(:,:,kk_sub,kk) end do ! end of kk_sub loop !----------------------------------------------------------------------- ! ! impose the boundary conditions by setting KAPPA=0 ! in the quarter cells adjacent to rigid boundaries. ! bottom B.C.s are considered within the kk-loop. ! !----------------------------------------------------------------------- ! B.C. at the bottom where (kk == KMT(:,:)) KAPPA_ISOP(:,:,kbt,kk) = 0.d0 KAPPA_THIC(:,:,kbt,kk) = 0.d0 end where enddo ! end of kk-loop ! B.C. at the top KAPPA_ISOP(:,:,ktp,1) = 0.d0 KAPPA_THIC(:,:,ktp,1) = 0.d0 FZTOP(:,:,:) = 0.d0 ! zero flux B.C. at the surface if ( transition_layer_on ) then call merged_streamfunction call apply_vertical_profile_to_isop_hor_diff endif endif ! end of k==1 if statement ! KMASK = merge(1.0d0, 0.d0, k < KMT(:,:)) KMASK = iw(:,:,k+1) !----------------------------------------------------------------------- ! ! calculate effective vertical diffusion coefficient ! NOTE: it is assumed that VDC has been set before this ! in vmix_coeffs or something similar. ! ! Dz(VDC * Dz(T)) where D is derivative rather than difference ! VDC = (Az(dz*Ax(KAPPA*HYX*SLX**2)) + Az(dz*Ay(KAPPA*HXY*SLY**2)))* ! dzw/TAREA ! !----------------------------------------------------------------------- if ( k < k1 ) then WORK1 = dzw(k)*KMASK*TAREA_R(:,:)* & ((1.0d0/odz(k))*0.25d0*KAPPA_ISOP(:,:,kbt,k)* & (HYX (:,:)*SLX(:,:,ieast, kbt,k)**2 & + HYXW(:,:)*SLX(:,:,iwest, kbt,k)**2 & + HXY (:,:)*SLY(:,:,jnorth,kbt,k)**2 & + HXYS(:,:)*SLY(:,:,jsouth,kbt,k)**2) & +(1.0d0/odz(k+1))*0.25d0*KAPPA_ISOP(:,:,ktp,k+1)* & (HYX (:,:)*SLX(:,:,ieast, ktp,k+1)**2 & + HYXW(:,:)*SLX(:,:,iwest, ktp,k+1)**2 & + HXY (:,:)*SLY(:,:,jnorth,ktp,k+1)**2 & + HXYS(:,:)*SLY(:,:,jsouth,ktp,k+1)**2)) do n=1,size(vdc,DIM=4) VDC_GM(:,:,k) = WORK1 vdc(:,:,k,n) = vdc(:,:,k,n) + WORK1 end do !! VDC_GM0(:,:,k) = merge(WORK1*dzwr(k), 0.d0, KMT(:,:) > k) ! VDC_GM0(:,:,k) = WORK1*dzwr(k)*iw(:,:,k+1) end if !----------------------------------------------------------------------- ! ! check if some horizontal diffusion needs to be added to the ! bottom half of the bottom cell ! !----------------------------------------------------------------------- if ( ah_bkg_bottom /= 0.d0 ) then where ( k == KMT(:,:) ) HOR_DIFF(:,:,kbt,k) = ah_bkg_bottom endwhere endif !----------------------------------------------------------------------- ! ! combine isopycnal and horizontal diffusion coefficients ! !----------------------------------------------------------------------- do j=1,j0 do i=1,i1 WORK3(i,j) = KAPPA_ISOP(i, j,ktp,k) & + HOR_DIFF (i, j,ktp,k) & + KAPPA_ISOP(i, j,kbt,k) & + HOR_DIFF (i, j,kbt,k) & + KAPPA_ISOP(i+1,j,ktp,k) & + HOR_DIFF (i+1,j,ktp,k) & + KAPPA_ISOP(i+1,j,kbt,k) & + HOR_DIFF (i+1,j,kbt,k) HOR_DIFF_X(i,j) = 0.25d0 * & ( HOR_DIFF (i, j,ktp,k) & + HOR_DIFF (i, j,kbt,k) & + HOR_DIFF (i+1,j,ktp,k) & + HOR_DIFF (i+1,j,kbt,k) ) KH_GM_X(i,j) = 0.25d0 * WORK3(i,j) enddo enddo do j=1,j1 do i=1,i0 WORK4(i,j) = KAPPA_ISOP(i,j, ktp,k) & + HOR_DIFF (i,j, ktp,k) & + KAPPA_ISOP(i,j, kbt,k) & + HOR_DIFF (i,j, kbt,k) & + KAPPA_ISOP(i,j+1,ktp,k) & + HOR_DIFF (i,j+1,ktp,k) & + KAPPA_ISOP(i,j+1,kbt,k) & + HOR_DIFF (i,j+1,kbt,k) HOR_DIFF_Y(i,j) = 0.25d0 * & ( HOR_DIFF (i,j, ktp,k) & + HOR_DIFF (i,j, kbt,k) & + HOR_DIFF (i,j+1,ktp,k) & + HOR_DIFF (i,j+1,kbt,k) ) KH_GM_Y(i,j) = 0.25d0 * WORK4(i,j) enddo enddo !----------------------------------------------------------------------- ! ! start loop over tracers ! !----------------------------------------------------------------------- kp1 = k + 1 if ( k == k1 ) kp1 = k if ( k < k1 ) then dz_bottom = 1.0d0/odz(kp1) factor = 1.0d0 else dz_bottom = 0.d0 factor = 0.d0 endif do n = 1,tnt !----------------------------------------------------------------------- ! ! calculate horizontal fluxes thru vertical faces of T-cell ! FX = dz*HYX*Ax(Az(KAPPA))*Dx(T) : flux in x-direction ! FY = dz*HXY*Ay(Az(KAPPA))*Dy(T) : flux in y-direction ! !----------------------------------------------------------------------- FX(:,:,n) = 1.0d0/odz(k) * CX * TDX(:,:,k,n) * WORK3 FY(:,:,n) = 1.0d0/odz(k) * CY * TDY(:,:,k,n) * WORK4 end do if ( .not. cancellation_occurs ) then !if(myid==0) write(*,*) .not. cancellation_occurs do j=1,j0 do i=1,i1 WORK1(i,j) = KAPPA_ISOP(i,j,ktp,k) & * SLX(i,j,ieast,ktp,k) * 1.0d0/odz(k) & - SF_SLX(i,j,ieast,ktp,k) WORK2(i,j) = KAPPA_ISOP(i,j,kbt,k) & * SLX(i,j,ieast,kbt,k) * 1.0d0/odz(k) & - SF_SLX(i,j,ieast,kbt,k) WORK3(i,j) = KAPPA_ISOP(i+1,j,ktp,k) & * SLX(i+1,j,iwest,ktp,k) * 1.0d0/odz(k) & - SF_SLX(i+1,j,iwest,ktp,k) WORK4(i,j) = KAPPA_ISOP(i+1,j,kbt,k) & * SLX(i+1,j,iwest,kbt,k) * 1.0d0/odz(k) & - SF_SLX(i+1,j,iwest,kbt,k) enddo enddo do n = 1,tnt if (n > 2 .and. k < k1) & TDZ(:,:,k+1,n) = TMIX(:,:,k ,n) - TMIX(:,:,k+1,n) do j=1,j0 do i=1,i1 FX0(i,j,n) = - CX(i,j) & * ( WORK1(i,j) * TDZ(i,j,k,n) * iw(i,j,k) & + WORK2(i,j) * TDZ(i,j,kp1,n) * iw(i,j,kp1) & + WORK3(i,j) * TDZ(i+1,j,k,n) * iw(i+1,j,k) & + WORK4(i,j) * TDZ(i+1,j,kp1,n) * iw(i+1,j,kp1) ) FX(i,j,n) = FX(i,j,n) + FX0(i,j,n) enddo enddo end do do j=1,j1 do i=1,i0 WORK1(i,j) = KAPPA_ISOP(i,j,ktp,k) & * SLY(i,j,jnorth,ktp,k) * 1.0d0/odz(k) & - SF_SLY(i,j,jnorth,ktp,k) WORK2(i,j) = KAPPA_ISOP(i,j,kbt,k) & * SLY(i,j,jnorth,kbt,k) * 1.0d0/odz(k) & - SF_SLY(i,j,jnorth,kbt,k) WORK3(i,j) = KAPPA_ISOP(i,j+1,ktp,k) & * SLY(i,j+1,jsouth,ktp,k) * 1.0d0/odz(k) & - SF_SLY(i,j+1,jsouth,ktp,k) WORK4(i,j) = KAPPA_ISOP(i,j+1,kbt,k) & * SLY(i,j+1,jsouth,kbt,k) * 1.0d0/odz(k) & - SF_SLY(i,j+1,jsouth,kbt,k) enddo enddo do n = 1,tnt do j=1,j1 do i=1,i0 FY0(i,j,n) = - CY(i,j) & * ( WORK1(i,j) * TDZ(i,j,k,n) * iw(i,j,k) & + WORK2(i,j) * TDZ(i,j,kp1,n) * iw(i,j,kp1) & + WORK3(i,j) * TDZ(i,j+1,k,n) * iw(i,j+1,k) & + WORK4(i,j) * TDZ(i,j+1,kp1,n) * iw(i,j+1,kp1) ) FY(i,j,n) = FY(i,j,n) + FY0(i,j,n) enddo enddo end do endif do n = 1,tnt !----------------------------------------------------------------------- ! ! calculate vertical fluxes thru horizontal faces of T-cell ! - Az(dz*Ax(HYX*KAPPA*SLX*TX)) - Az(dz*Ay(HXY*KAPPA*SLY*TY)) ! calculate isopycnal diffusion from flux differences ! DTK = (Dx(FX)+Dy(FY)+Dz(FZ)) / volume ! !----------------------------------------------------------------------- ! GTK(:,:,n) = 0.d0 GTK0(:,:,n) = 0.d0 if ( k < k1 ) then WORK3 = 0.d0 if ( .not. cancellation_occurs ) then !if(myid==0) write(*,*) .not. cancellation_occurs !pw loop split to improve performance -- 2 do j=2,j1 do i=2,i1 WORK3(i,j) = WORK3(i,j) & + ( (1.0d0/odz(k)) * KAPPA_ISOP(i,j,kbt,k) & * ( SLX(i,j,ieast ,kbt,k) & * HYX(i ,j) * TDX(i ,j ,k ,n) & + SLY(i,j,jnorth,kbt,k) & * HXY(i ,j) * TDY(i ,j ,k ,n) & + SLX(i,j,iwest ,kbt,k) & * HYX(i-1,j) * TDX(i-1,j ,k ,n) & + SLY(i,j,jsouth,kbt,k) & * HXY(i ,j-1) * TDY(i ,j-1,k ,n) ) ) enddo enddo do j=2,j1 do i=2,i1 WORK3(i,j) = WORK3(i,j) & + ( SF_SLX(i ,j ,ieast ,kbt,k) & * HYX(i ,j) * TDX(i ,j ,k ,n) & + SF_SLY(i ,j ,jnorth,kbt,k) & * HXY(i ,j) * TDY(i ,j ,k ,n) & + SF_SLX(i ,j ,iwest ,kbt,k) & * HYX(i-1,j) * TDX(i-1,j ,k ,n) & + SF_SLY(i ,j ,jsouth,kbt,k) & * HXY(i ,j-1) * TDY(i ,j-1,k ,n) ) enddo enddo do j=2,j1 do i=2,i1 WORK3(i,j) = WORK3(i,j) & + ( dz_bottom * KAPPA_ISOP(i,j,ktp,kp1) & * ( SLX(i ,j ,ieast ,ktp,kp1) & * HYX(i ,j) * TDX(i ,j ,kp1,n) & + SLY(i ,j ,jnorth,ktp,kp1) & * HXY(i ,j) * TDY(i ,j ,kp1,n) & + SLX(i ,j ,iwest ,ktp,kp1) & * HYX(i-1,j) * TDX(i-1,j ,kp1,n) & + SLY(i ,j ,jsouth,ktp,kp1) & * HXY(i ,j-1) * TDY(i ,j-1,kp1,n) ) ) enddo enddo do j=2,j1 do i=2,i1 WORK3(i,j) = WORK3(i,j) & + ( factor & * ( SF_SLX(i ,j ,ieast ,ktp,kp1) & * HYX(i ,j) * TDX(i ,j ,kp1,n) & + SF_SLY(i ,j ,jnorth,ktp,kp1) & * HXY(i ,j) * TDY(i ,j ,kp1,n) & + SF_SLX(i ,j ,iwest ,ktp,kp1) & * HYX(i-1,j) * TDX(i-1,j ,kp1,n) & + SF_SLY(i ,j ,jsouth,ktp,kp1) & * HXY(i ,j-1) * TDY(i ,j-1,kp1,n) ) ) enddo enddo !if(myid.eq.0)write(*,103)'myid0_x_FZTOP',sum(FZTOP(2:i1,2:j1,1)),k!FX0(42,37,1),FX0(i1,37,1),FX0(i0,37,1),F !103 FORMAT(1X,F15.6) do j=2,j1 do i=2,i1 fz = -KMASK(i,j) * 0.25d0 * WORK3(i,j) ! GTK(i,j,n) = ( FX(i,j,n) - FX(i-1,j,n) & ! + FY(i,j,n) - FY(i,j-1,n) & ! + FZTOP(i,j,n) - fz )*odz(k)*TAREA_R(i,j) GTK0(i,j,n) = (FX0(i,j,n) - FX0(i-1,j,n) & + FY0(i,j,n) - FY0(i,j-1,n) & + FZTOP(i,j,n) - fz )*odz(k)*TAREA_R(i,j) !!! GTK0(i,j,n) = FX0(i,j,n)*odz(k)*TAREA_R(i,j) -FX0(i-1,j,n)*odz(k)*TAREA_R(i,j) & !!! + FY0(i,j,n)*odz(k)*TAREA_R(i,j) -FY0(i,j-1,n)*odz(k)*TAREA_R(i,j) & !!! +FZTOP(i,j,n)*odz(k)*TAREA_R(i,j) -fz*odz(k)*TAREA_R(i,j) FZTOP(i,j,n) = fz enddo enddo !if (myid.eq.0)write(*,101)'myid0_x_fz',sum(FZTOP(2:i1,2:j1,1)),k!FX0(42,37,1),FX0(i1,37,1),FX0(i0,37,1),FX0(i1,j1,1) !101 FORMAT(1X,F15.6) !if(myid.eq.0)write(*,111)'myid0_x_fz2',FZTOP(2,2,1),k!FX0(42,37,1),FX0(i1,37,1),FX0(i0,37,1),FX0(i1,j1,1) !111 FORMAT(1X,F15.6) !if (myid.eq.8) !write(*,*)'myid8_x_KMASK',FZTOP(2,1,1),FZTOP(2,37,1),k!FX0(42,37,1),FX0(i1,37,1),FX0(i0,37,1),FX0(i1,j1,1)G !GGG(k)=sum(GTK0(2:i1,2:j1,1)) !if (myid.eq.0)write(*,102)'myid0_x_GTK0',sum(GTK0(2:i1,2:j1,1)),k!FX0(42,37,1),FX0(i1,37,1),FX0(i0,37,1),FX0(i1,j1 !102 FORMAT(1X,F15.6) else !pw loop split to improve performance do j=2,j1 do i=2,i1 WORK3(i,j) = & ( (1.0d0/odz(k)) * KAPPA_ISOP(i,j,kbt,k) & * ( SLX(i,j,ieast ,kbt,k) & * HYX(i ,j) * TDX(i ,j ,k ,n) & + SLY(i,j,jnorth,kbt,k) & * HXY(i ,j) * TDY(i ,j ,k ,n) & + SLX(i,j,iwest ,kbt,k) & * HYX(i-1,j) * TDX(i-1,j ,k ,n) & + SLY(i,j,jsouth,kbt,k) & * HXY(i ,j-1) * TDY(i ,j-1,k ,n) ) ) enddo enddo do j=2,j1 do i=2,i1 WORK3(i,j) = WORK3(i,j) & + ( dz_bottom * KAPPA_ISOP(i,j,ktp,kp1) & * ( SLX(i ,j ,ieast ,ktp,kp1) & * HYX(i ,j) * TDX(i ,j ,kp1,n) & + SLY(i ,j ,jnorth,ktp,kp1) & * HXY(i ,j) * TDY(i ,j ,kp1,n) & + SLX(i ,j ,iwest ,ktp,kp1) & * HYX(i-1,j) * TDX(i-1,j ,kp1,n) & + SLY(i ,j ,jsouth,ktp,kp1) & * HXY(i ,j-1) * TDY(i ,j-1,kp1,n) ) ) enddo enddo do j=2,j1 do i=2,i1 fz = -KMASK(i,j) * 0.25d0 * WORK3(i,j) ! GTK(i,j,n) = ( FX(i,j,n) - FX(i-1,j,n) & ! + FY(i,j,n) - FY(i,j-1,n) & ! + FZTOP(i,j,n) - fz )*odz(k)*TAREA_R(i,j) GTK0(i,j,n) = ( FZTOP(i,j,n) - fz )*odz(k)*TAREA_R(i,j) FZTOP(i,j,n) = fz enddo enddo endif else ! k = km !if (myid.eq.0) write(*,100)'myid0_x_FZTOPbottom',sum(FZTOP(:,:,1)),k!FX0(42,37,1),FX0(i1,37,1),FX0(i0,37,1),F !100 FORMAT(1X,F15.6) do j=2,j1 do i=2,i1 ! GTK(i,j,n) = ( FX(i,j,n) - FX(i-1,j,n) & ! + FY(i,j,n) - FY(i,j-1,n) & ! + FZTOP(i,j,n) )*odz(k)*TAREA_R(i,j) if ( .not. cancellation_occurs ) then GTK0(i,j,n) = (FX0(i,j,n) - FX0(i-1,j,n) & + FY0(i,j,n) - FY0(i,j-1,n) & + FZTOP(i,j,n) )*odz(k)*TAREA_R(i,j) !!! GTK0(i,j,n) = FX0(i,j,n)*odz(k)*TAREA_R(i,j)-FX0(i-1,j,n)*odz(k)*TAREA_R(i,j) & !!! + FY0(i,j,n)*odz(k)*TAREA_R(i,j)-FY0(i,j-1,n)*odz(k)*TAREA_R(i,j) & !!! + FZTOP(i,j,n)*odz(k)*TAREA_R(i,j) else GTK0(i,j,n) = FZTOP(i,j,n)*odz(k)*TAREA_R(i,j) endif FZTOP(i,j,n) = 0.d0 enddo enddo endif !----------------------------------------------------------------------- ! ! end of tracer loop ! !----------------------------------------------------------------------- !TS ! GTK0(2:i1,2:j1,n)=GTK0(2:i1,2:j1,n)*IN(2:i1,2:j1,k)*& ! IN(1:i2,2:j1,k)*IN(3:i0,2:j1,k)*IN(2:i1,1:j2,k)*IN(2:i1,3:j0,k)*& ! IN(1:i2,2:j1,kp1)*IN(3:i0,2:j1,kp1)*IN(2:i1,1:j2,kp1)*IN(2:i1,3:j0,kp1) enddo !----------------------------------------------------------------------- ! ! diagnostic computation of the bolus velocities ! !----------------------------------------------------------------------- if ( diag_gm_bolus ) then do j=1,j1 do i=1,i1 WORK1(i,j) = ( SF_SLX(i ,j,ieast,kbt,k) & + factor * SF_SLX(i ,j,ieast,ktp,kp1) & + SF_SLX(i+1,j,iwest,kbt,k) & + factor * SF_SLX(i+1,j,iwest,ktp,kp1)) & * 0.25d0 * HYX(i,j) WORK2(i,j) = ( SF_SLY(i,j ,jnorth,kbt,k) & + factor * SF_SLY(i,j ,jnorth,ktp,kp1) & + SF_SLY(i,j+1,jsouth,kbt,k) & + factor * SF_SLY(i,j+1,jsouth,ktp,kp1)) & * 0.25d0 * HXY(i,j) enddo enddo ! UIB = merge( WORK1, 0.d0, k < KMT(:,:) & ! .and. k < KMTE(:,:) ) ! VIB = merge( WORK2, 0.d0, k < KMT(:,:) & ! .and. k < KMTN(:,:) ) if ( k < k1 ) then UIB = WORK1*iu(1:i0,:,k+1) VIB = WORK2*iv(:,1:j0,k+1) else UIB = 0.d0 VIB = 0.d0 endif ! WORK1 = merge( UIT(:,:) - UIB, 0.d0, k <= KMT(:,:) & ! .and. k <= KMTE(:,:) ) ! WORK2 = merge( VIT(:,:) - VIB, 0.d0, k <= KMT(:,:) & ! .and. k <= KMTN(:,:) ) if (k==1) then WORK1 = (UIT-UIB)*iu0(1:i0,:) WORK2 = (VIT-VIB)*iv0(:,1:j0) else WORK1 = (UIT-UIB)*iu(1:i0,:,k) WORK2 = (VIT-VIB)*iv(:,1:j0,k) endif U_ISOP = WORK1 * odz(k) / HTE V_ISOP = WORK2 * odz(k) / HTN if ( linertial .and. k == 2 ) then bolus_sp(:,:) = 50.0d0* sqrt(U_ISOP**2.0d0 + V_ISOP**2.0d0) endif do j=2,j1 do i=2,i1 if ( k < KMT(i,j) ) then WBOT_ISOP(i,j) = WTOP_ISOP(i,j)+ TAREA_R(i,j) & * ( WORK1(i,j) - WORK1(i-1,j) & + WORK2(i,j) - WORK2(i,j-1) ) else WBOT_ISOP(i,j) = 0.d0 endif enddo enddo endif !----------------------------------------------------------------------- ! ! update remaining bottom-face fields to top-face fields for next ! pass ! !----------------------------------------------------------------------- if ( diag_gm_bolus ) then UIT(:,:) = UIB VIT(:,:) = VIB endif !----------------------------------------------------------------------- ! ! accumulate time average if necessary; testing is internal to ! accumulate_tavg_field ! !----------------------------------------------------------------------- if ( diag_gm_bolus ) then UISOP(:,:,k) = U_ISOP VISOP(:,:,k) = V_ISOP WISOP(:,:,k) = WTOP_ISOP if (ldiag_budgets) then ! tavg_ADVT_ISOP ! Vertically-Integrated T Eddy-Induced Advection Tendency (diagnostic) ! ADVT_ISOP [UNIT:cm degC/s] WORK1 = 0.5d0 * HTE(:,:) * U_ISOP * ( TMIX(:,:,k,1) & + eoshift(TMIX(:,:,k,1), dim=1, shift=1) ) WORK2 = eoshift(WORK1, dim=1, shift=-1) WORK3 = WORK1 - WORK2 WORK1 = 0.5d0 * HTN(:,:) * V_ISOP * ( TMIX(:,:,k,1) & + eoshift(TMIX(:,:,k,1), dim=2, shift=1) ) WORK2 = eoshift(WORK1, dim=2, shift=-1) WORK3 = WORK3 + WORK1 - WORK2 WORK1 = 0.0d0 if (k==1) ADVT_ISOP = 0.0d0 do j=2,j1 do i=2,i1 if ( k <= KMT(i,j) ) then WORK1(i,j) = - (1.0d0/odz(k)) * TAREA_R(i,j) * WORK3(i,j) endif enddo enddo ADVT_ISOP = ADVT_ISOP + WORK1 ! tavg_ADVS_ISOP ! Vertically-Integrated S Eddy-Induced Advection Tendency (diagnostic) ! ADVS_ISOP [UNIT:cm gram/kilogram/s] WORK1 = 0.50d0 * HTE(:,:) * U_ISOP * ( TMIX(:,:,k,2) & + eoshift(TMIX(:,:,k,2), dim=1, shift=1) ) WORK2 = eoshift(WORK1, dim=1, shift=-1) WORK3 = WORK1 - WORK2 WORK1 = 0.50d0 * HTN(:,:) * V_ISOP * ( TMIX(:,:,k,2) & + eoshift(TMIX(:,:,k,2), dim=2, shift=1) ) WORK2 = eoshift(WORK1, dim=2, shift=-1) WORK3 = WORK3 + WORK1 - WORK2 WORK1 = 0.0d0 if (k==1) ADVS_ISOP = 0.0d0 do j=2,j1 do i=2,i1 if ( k <= KMT(i,j) ) then WORK1(i,j) = - (1.0d0/odz(k)) * TAREA_R(i,j) * WORK3(i,j) & *1000.0d0 endif enddo enddo ADVS_ISOP = ADVS_ISOP + WORK1 ! tavg_VNT_ISOP ! Heat Flux Tendency in grid-y Dir due to Eddy-Induced Vel (diagnostic) ! VNT_ISOP [UNIT:degC/s] WORK1 = 0.5d0 * V_ISOP * HTN(:,:) * TAREA_R(:,:) VNT_ISOP(:,:,k) = WORK1 * ( TMIX(:,:,k,1) & + eoshift(TMIX(:,:,k,1), dim=2, shift=1) ) ! tavg_VNS_ISOP ! Salt Flux Tendency in grid-y Dir due to Eddy-Induced Vel (diagnostic) ! VNS_ISOP [UNIT:gram/kilogram/s] VNS_ISOP(:,:,k) = WORK1 * ( TMIX(:,:,k,2) & + eoshift(TMIX(:,:,k,2), dim=2, shift=1) )*1000.0d0 endif ! ldiag_budgets on endif ! bolus velocity option on if (ldiag_budgets) then do n = 1,tnt ! tavg_HDIFE_TRACER(n) ! Horizontal Diffusive Flux in grid-x direction ! HDIFE_TRACER(n)[TEMPdegC/s][SALTgram/kilogram/s] do j=2,j1 do i=2,i1 HDIFE(i,j,k,n) = FX(i,j,n)*odz(k)*TAREA_R(i,j) if (n==2) HDIFE(i,j,k,n) = HDIFE(i,j,k,n) *1000.0d0 enddo enddo ! tavg_HDIFN_TRACER(n) ! Horizontal Diffusive Flux in grid-y direction ! HDIFN_TRACER(n)[TEMPdegC/s][SALTgram/kilogram/s] do j=2,j1 do i=2,i1 HDIFN(i,j,k,n) = FY(i,j,n)*odz(k)*TAREA_R(i,j) if (n==2) HDIFN(i,j,k,n) = HDIFN(i,j,k,n)*1000.0d0 enddo enddo ! tavg_HDIFB_TRACER(n) ! Horizontal Diffusive Flux across Bottom Face ! tavg_HDIFB_TRACER(n)[TEMPdegC/s][SALTgram/kilogram/s] do j=2,j1 do i=2,i1 HDIFB(i,j,k,n) = FZTOP(i,j,n)*odz(k)*TAREA_R(i,j) if (n==2) HDIFB(i,j,k,n) = HDIFB(i,j,k,n)*1000.0d0 enddo enddo enddo endif ! ldiag_budgets on GTK0_T(:,:,k) = GTK0(:,:,1) !UNIT=degree C/[T] GTK0_S(:,:,k) = GTK0(:,:,2)*1000.0d0 !UNIT=(msu/[T]=g/g/[T])*1000 ! = psu/[T]=g/kg/[T] !----------------------------------------------------------------------- !EOC end subroutine hdifft_gm !*********************************************************************** !BOP ! !IROUTINE: buoyancy_frequency_dependent_profile ! !INTERFACE: subroutine buoyancy_frequency_dependent_profile (TMIX) ! !DESCRIPTION: ! Computes normalized buoyancy frequency ($N$) dependent profiles that ! are used to represent the vertical variation of the isopycnal and/or ! thickness diffusion coefficients. We use ! \begin{equation} ! KAPPA_VERTICAL = {{N^2}\over {N_ref^2}}, ! \end{equation} ! where $N_ref$ is the value of $N$ "just" below the surface diabatic ! layer (SDL) provided that $N^2 > 0$ there. Otherwise, $N_ref$ is the ! first stable $N^2$ below SDL. KAPPA_VERTICAL is set to unity for shallower ! depths. Also, 0.1 <= KAPPA_VERTICAL <= 1.0 is enforced. Note that ! this implies KAPPA_VERTICAL = 0.1 if the local $N^2 \le 0$. ! KAPPA_VERTICAL is located at the model tracer points. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: real*8, dimension(i0,j0,k1,tnt), intent(in) :: & TMIX ! tracers at all vertical levels ! at mixing time level !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer :: k ! vertical loop index integer, dimension(i0,j0) :: & K_MIN ! k index below SDL real*8, dimension(i0,j0) :: & TEMP_K, & ! temperature at level k TEMP_KP1, & ! temperature at level k+1 RHOT, & ! dRHO/dT RHOS, & ! dRHO/dS BUOY_FREQ_SQ_REF, & ! reference (normalization) value ! of N^2 SDL ! surface diabatic layer (see below) real*8, dimension(i0,j0,k1) :: & BUOY_FREQ_SQ_NORM ! normalized N^2 defined at level interfaces !----------------------------------------------------------------------- ! ! initialize variables ! ! note that "km+1" value in K_MIN is used as a flag, i.e. it indicates ! that a minimum value has not been found yet or that it does not ! exist. ! !----------------------------------------------------------------------- BUOY_FREQ_SQ_NORM = 0.d0 BUOY_FREQ_SQ_REF = 0.d0 KAPPA_VERTICAL(:,:,:) = 1.0d0 K_MIN = merge( k1+1, 0, KMT(:,:) /= 0 ) SDL = kpp_hblt(:,:) if ( transition_layer_on ) SDL = TLT%INTERIOR_DEPTH(:,:) !----------------------------------------------------------------------- ! ! compute N^2 ! !----------------------------------------------------------------------- do k=1,k1-1 if ( k == 1 ) & TEMP_K = max( -2.0d0, TMIX(:,:,k,1) ) TEMP_KP1 = max( -2.0d0, TMIX(:,:,k+1,1) ) call state( k, k+1, TMIX(:,:,k,1), TMIX(:,:,k,2), & DRHODT=RHOT, DRHODS=RHOS ) where ( k < KMT(:,:) ) BUOY_FREQ_SQ(:,:,k) = max( 0.d0, - grav * dzwr(k) & * ( RHOT * ( TEMP_K - TEMP_KP1 ) & + RHOS * ( TMIX(:,:,k, 2) - TMIX(:,:,k+1,2) ) ) ) endwhere TEMP_K = TEMP_KP1 enddo !----------------------------------------------------------------------- ! ! determine the reference buoyancy frequency and the associated ! level index (most likely just below SDL) ! !----------------------------------------------------------------------- do k=1,k1-1 where ( ( K_MIN == k1+1 ) .and. ( (z(2*k+1) - z(1)) > SDL ) .and. & ( k <= KMT(:,:) ) .and. & ( BUOY_FREQ_SQ(:,:,k) > 0.d0 ) ) BUOY_FREQ_SQ_REF = BUOY_FREQ_SQ(:,:,k) K_MIN = k endwhere enddo !----------------------------------------------------------------------- ! ! now compute the normalized profiles at the interfaces ! !----------------------------------------------------------------------- do k=1,k1-1 where ( ( k >= K_MIN ) .and. ( k < KMT(:,:) ) .and. & ( BUOY_FREQ_SQ_REF /= 0.d0 ) ) BUOY_FREQ_SQ_NORM(:,:,k) = & max( BUOY_FREQ_SQ(:,:,k) / BUOY_FREQ_SQ_REF, 0.1d0 ) BUOY_FREQ_SQ_NORM(:,:,k) = & min( BUOY_FREQ_SQ_NORM(:,:,k), 1.0d0 ) elsewhere BUOY_FREQ_SQ_NORM(:,:,k) = 1.0d0 endwhere enddo do k=1,k1-1 where ( k == KMT(:,:)-1 ) BUOY_FREQ_SQ_NORM(:,:,k+1) = BUOY_FREQ_SQ_NORM(:,:,k) endwhere enddo !----------------------------------------------------------------------- ! ! finally, do not average interface values of BUOY_FREQ_SQ_NORM to ! the tracer points, but instead copy from above to preserve the ! extrema. ! !----------------------------------------------------------------------- do k=2,k1 where ( ( k > K_MIN ) .and. ( k <= KMT(:,:) ) ) KAPPA_VERTICAL(:,:,k) = BUOY_FREQ_SQ_NORM(:,:,k-1) endwhere enddo !----------------------------------------------------------------------- !EOC end subroutine buoyancy_frequency_dependent_profile !*********************************************************************** !BOP ! !IROUTINE: transition_layer ! !INTERFACE: subroutine transition_layer ! !DESCRIPTION: ! Compute transition layer related fields. the related algorithms ! should work even with zero transition layer depth. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer :: & k, kk ! loop indices integer, dimension(i0,j0) :: & K_START, & ! work arrays for TLT%K_LEVEL and K_SUB ! TLT%ZTW, respectively logical, dimension(i0,j0) :: & COMPUTE_TLT ! flag real*8, dimension(i0,j0) :: & WORK ! work space for TLT%THICKNESS real*8, dimension(2) :: & reference_depth ! zt or zw !----------------------------------------------------------------------- ! ! initialize various quantities ! !----------------------------------------------------------------------- K_START = 0 K_SUB = 0 TLT%THICKNESS(:,:) = 0.d0 TLT%INTERIOR_DEPTH(:,:) = 0.d0 TLT%K_LEVEL(:,:) = 0 TLT%ZTW(:,:) = 0 COMPUTE_TLT = merge(.true., .false., KMT(:,:) /= 0) !----------------------------------------------------------------------- ! ! initial pass to determine the minimum transition layer thickness. ! the vertical level at or below which the interior region starts ! is also computed. ! !----------------------------------------------------------------------- do k=1,k1 where ( COMPUTE_TLT .and. & TLT%DIABATIC_DEPTH(:,:) < (z(2*k+1) - z(1)) ) K_START = k+1 K_SUB = ktp TLT%THICKNESS(:,:) = (z(2*k+1) - z(1)) - TLT%DIABATIC_DEPTH(:,:) TLT%K_LEVEL(:,:) = k TLT%ZTW(:,:) = 2 COMPUTE_TLT = .false. endwhere where ( k /= 1 .and. K_START == k+1 .and. & TLT%DIABATIC_DEPTH(:,:) < (z(2*k) - z(1)) ) K_START = k K_SUB = kbt TLT%THICKNESS(:,:) = (z(2*k) - z(1)) - TLT%DIABATIC_DEPTH(:,:) TLT%K_LEVEL(:,:) = k TLT%ZTW(:,:) = 1 endwhere enddo !----------------------------------------------------------------------- ! ! using R * |S| as the vertical displacement length scale associated ! with a horizontal displacement equivalent to the Rossby deformation ! radius, R, determine if the transition layer thickness needs to be ! modified. |S| is the absolute value of the slope of the isopycnal ! surfaces. First, start with columns where K_SUB = kbt. ! !----------------------------------------------------------------------- where ( KMT(:,:) == 0 .or. K_START > KMT(:,:) .or. & ( K_START == KMT(:,:) .and. K_SUB == kbt ) ) COMPUTE_TLT = .false. elsewhere COMPUTE_TLT = .true. endwhere do k=1,k1-1 WORK = 0.d0 where ( COMPUTE_TLT .and. K_SUB == kbt .and. & K_START < KMT(:,:) .and. K_START == k ) WORK = max(SLA_SAVE(:,:,kbt,k), & SLA_SAVE(:,:,ktp,k+1)) * RB(:,:) endwhere where ( WORK /= 0.d0 .and. & TLT%DIABATIC_DEPTH(:,:) < ((z(2*k+1) - z(1)) - WORK) ) COMPUTE_TLT = .false. endwhere where ( WORK /= 0.d0 .and. & TLT%DIABATIC_DEPTH(:,:) >= ((z(2*k+1) - z(1)) - WORK) ) K_START = K_START + 1 K_SUB = ktp TLT%THICKNESS(:,:) = (z(2*k+1) - z(1)) - TLT%DIABATIC_DEPTH(:,:) TLT%K_LEVEL(:,:) = k TLT%ZTW(:,:) = 2 endwhere enddo !----------------------------------------------------------------------- ! ! now consider the deeper levels ! !----------------------------------------------------------------------- do k=2,k1 reference_depth(ktp) = z(2*k) - z(1) reference_depth(kbt) = z(2*k+1) - z(1) do kk=ktp,kbt WORK = 0.d0 if (kk == ktp) then where ( COMPUTE_TLT .and. K_START <= KMT(:,:) .and. & K_START == k ) WORK = max(SLA_SAVE(:,:,ktp,k), & SLA_SAVE(:,:,kbt,k)) * RB(:,:) endwhere else ! Checking k < km guarantees that k+1 is not out of bounds if (k.lt.k1) then where ( COMPUTE_TLT .and. K_START < KMT(:,:) .and. & K_START == k ) WORK = max(SLA_SAVE(:,:,kbt,k), & SLA_SAVE(:,:,ktp,k+1)) * RB(:,:) endwhere end if where ( COMPUTE_TLT .and. K_START == KMT(:,:) .and. & K_START == k ) WORK = SLA_SAVE(:,:,kbt,k) * RB(:,:) endwhere endif where ( WORK /= 0.d0 .and. & TLT%DIABATIC_DEPTH(:,:) < (reference_depth(kk) - WORK) ) COMPUTE_TLT = .false. endwhere where ( WORK /= 0.d0 .and. & TLT%DIABATIC_DEPTH(:,:) >= (reference_depth(kk) - WORK) ) TLT%THICKNESS(:,:) = reference_depth(kk) & - TLT%DIABATIC_DEPTH(:,:) TLT%K_LEVEL(:,:) = k TLT%ZTW(:,:) = kk endwhere enddo where ( COMPUTE_TLT .and. K_START == k ) K_START = K_START + 1 endwhere enddo !----------------------------------------------------------------------- ! ! compute the depth at which the interior, adiabatic region starts ! !----------------------------------------------------------------------- do k=1,k1 where ( TLT%K_LEVEL(:,:) == k .and. & TLT%ZTW(:,:) == 1 ) TLT%INTERIOR_DEPTH(:,:) = z(2*k) - z(1) endwhere where ( TLT%K_LEVEL(:,:) == k .and. & TLT%ZTW(:,:) == 2 ) TLT%INTERIOR_DEPTH(:,:) = z(2*k+1) - z(1) endwhere enddo COMPUTE_TLT = .false. where ( KMT(:,:) /= 0 .and. & TLT%INTERIOR_DEPTH(:,:) == 0.d0 ) & COMPUTE_TLT = .true. where ( KMT(:,:) == 0 .and. & TLT%INTERIOR_DEPTH(:,:) /= 0.d0 ) & COMPUTE_TLT = .true. !----------------------------------------------------------------------- !EOC end subroutine transition_layer !*********************************************************************** !BOP ! !IROUTINE: merged_streamfunction ! !INTERFACE: subroutine merged_streamfunction ! !DESCRIPTION: ! Construct a merged streamfunction that has the appropriate ! behavior in the surface diabatic region, transition layer, and ! adiabatic interior ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer :: k, kk ! loop indices real*8, dimension(i0,j0,2) :: & WORK1, WORK2, WORK3, WORK4 ! work arrays real*8, dimension(i0,j0) :: & WORK2_NEXT, WORK4_NEXT ! WORK2 or WORK4 at next level real*8, dimension(i0,j0) :: & WORK5, WORK6, WORK7 ! more work arrays logical, dimension(i0,j0) :: & LMASK ! flag real*8, dimension(2) :: & reference_depth ! zt or zw !----------------------------------------------------------------------- ! ! initialize various quantities ! !----------------------------------------------------------------------- SF_SLX(:,:,:,:,:) = 0.d0 SF_SLY(:,:,:,:,:) = 0.d0 WORK1 = 0.d0 WORK2 = 0.d0 WORK3 = 0.d0 WORK4 = 0.d0 WORK5 = 0.d0 WORK6 = 0.d0 WORK7 = 0.d0 WORK2_NEXT = 0.d0 WORK4_NEXT = 0.d0 !----------------------------------------------------------------------- ! ! compute the interior streamfunction and its first derivative at the ! INTERIOR_DEPTH level. WORK1 and WORK2 contain the streamfunction ! and its first derivative, respectively, for the zonal component ! for the east and west sides of a grid cell. WORK3 and WORK4 are ! the corresponding fields for the meridional component for the ! north and south sides of a grid cell. Note that these definitions ! include a "dz". Also, the first derivative computations assume ! that the streamfunctions are located in the middle of the top or ! bottom half of a grid cell, hence a factor of two in WORK2 and ! WORK4 calculations. ! !----------------------------------------------------------------------- do k=1,k1-1 do kk=1,2 LMASK = TLT%K_LEVEL(:,:) == k .and. & TLT%K_LEVEL(:,:) < KMT(:,:) .and. & TLT%ZTW(:,:) == 1 where ( LMASK ) WORK1(:,:,kk) = KAPPA_THIC(:,:,kbt,k) & * SLX(:,:,kk,kbt,k) * (1.0d0/odz(k)) WORK2(:,:,kk) = 2.0d0 * dzwr(k) * ( WORK1(:,:,kk) & - KAPPA_THIC(:,:,ktp,k+1) * SLX(:,:,kk,ktp,k+1) & * (1.0d0/odz(k+1)) ) WORK2_NEXT = 2.0d0 * ( & KAPPA_THIC(:,:,ktp,k+1) * SLX(:,:,kk,ktp,k+1) - & KAPPA_THIC(:,:,kbt,k+1) * SLX(:,:,kk,kbt,k+1) ) WORK3(:,:,kk) = KAPPA_THIC(:,:,kbt,k) & * SLY(:,:,kk,kbt,k) * (1.0d0/odz(k)) WORK4(:,:,kk) = 2.0d0 * dzwr(k) * ( WORK3(:,:,kk) & - KAPPA_THIC(:,:,ktp,k+1) * SLY(:,:,kk,ktp,k+1) & * (1.0d0/odz(k+1)) ) WORK4_NEXT = 2.0d0 * ( & KAPPA_THIC(:,:,ktp,k+1) * SLY(:,:,kk,ktp,k+1) - & KAPPA_THIC(:,:,kbt,k+1) * SLY(:,:,kk,kbt,k+1) ) endwhere where ( LMASK .and. abs(WORK2_NEXT) < abs(WORK2(:,:,kk)) ) & WORK2(:,:,kk) = WORK2_NEXT where ( LMASK .and. abs(WORK4_NEXT) < abs(WORK4(:,:,kk)) ) & WORK4(:,:,kk) = WORK4_NEXT LMASK = TLT%K_LEVEL(:,:) == k .and. & TLT%K_LEVEL(:,:) < KMT(:,:) .and. & TLT%ZTW(:,:) == 2 where ( LMASK ) WORK1(:,:,kk) = KAPPA_THIC(:,:,ktp,k+1) & * SLX(:,:,kk,ktp,k+1) WORK2(:,:,kk) = 2.0d0 * ( WORK1(:,:,kk) & - ( KAPPA_THIC(:,:,kbt,k+1) & * SLX(:,:,kk,kbt,k+1) ) ) WORK1(:,:,kk) = WORK1(:,:,kk) * (1.0d0/odz(k+1)) WORK3(:,:,kk) = KAPPA_THIC(:,:,ktp,k+1) & * SLY(:,:,kk,ktp,k+1) WORK4(:,:,kk) = 2.0d0 * ( WORK3(:,:,kk) & - ( KAPPA_THIC(:,:,kbt,k+1) & * SLY(:,:,kk,kbt,k+1) ) ) WORK3(:,:,kk) = WORK3(:,:,kk) * (1.0d0/odz(k+1)) endwhere LMASK = LMASK .and. TLT%K_LEVEL(:,:) + 1 < KMT(:,:) if (k.lt.k1-1) then ! added to avoid out of bounds access where ( LMASK ) WORK2_NEXT = 2.0d0 * dzwr(k+1) * ( & KAPPA_THIC(:,:,kbt,k+1) * SLX(:,:,kk,kbt,k+1) & * (1.0d0/odz(k+1))- & KAPPA_THIC(:,:,ktp,k+2) * SLX(:,:,kk,ktp,k+2) & * (1.0d0/odz(k+2))) WORK4_NEXT = 2.0d0 * dzwr(k+1) * ( & KAPPA_THIC(:,:,kbt,k+1) * SLY(:,:,kk,kbt,k+1) & * (1.0d0/odz(k+1))- & KAPPA_THIC(:,:,ktp,k+2) * SLY(:,:,kk,ktp,k+2) & * (1.0d0/odz(k+2))) endwhere end if where ( LMASK .and. abs(WORK2_NEXT) < abs(WORK2(:,:,kk)) ) & WORK2(:,:,kk) = WORK2_NEXT where ( LMASK .and. abs(WORK4_NEXT) < abs(WORK4(:,:,kk)) ) & WORK4(:,:,kk) = WORK4_NEXT enddo enddo !----------------------------------------------------------------------- ! ! compute the depth independent interpolation factors used in the ! linear and quadratic interpolations within the diabatic and ! transition regions, respectively. ! !----------------------------------------------------------------------- WORK5 = 0.d0 where (KMT(:,:) /= 0) WORK5(:,:) = 1.0d0 / ( 2.0d0 * TLT%DIABATIC_DEPTH(:,:) & + TLT%THICKNESS(:,:) ) endwhere WORK6 = 0.d0 where ((KMT(:,:) /= 0) .AND. (TLT%THICKNESS(:,:) > eps)) WORK6(:,:) = WORK5(:,:) / TLT%THICKNESS(:,:) endwhere !----------------------------------------------------------------------- ! ! start of interpolation to construct the merged streamfunction ! !----------------------------------------------------------------------- do k=1,k1 reference_depth(ktp) = (z(2*k) - z(1)) - 0.25d0 * (1.0d0/odz(k)) reference_depth(kbt) = (z(2*k) - z(1)) + 0.25d0 * (1.0d0/odz(k)) do kk=ktp,kbt !----------------------------------------------------------------------- ! ! diabatic region: use linear interpolation (in streamfunction) ! !----------------------------------------------------------------------- where ( reference_depth(kk) <= TLT%DIABATIC_DEPTH(:,:) & .and. k <= KMT(:,:) ) SF_SLX(:,:,1,kk,k) = reference_depth(kk) * WORK5 & * ( 2.0d0 * WORK1(:,:,1) + TLT%THICKNESS(:,:) & * WORK2(:,:,1) ) SF_SLX(:,:,2,kk,k) = reference_depth(kk) * WORK5 & * ( 2.0d0 * WORK1(:,:,2) + TLT%THICKNESS(:,:) & * WORK2(:,:,2) ) SF_SLY(:,:,1,kk,k) = reference_depth(kk) * WORK5 & * ( 2.0d0 * WORK3(:,:,1) + TLT%THICKNESS(:,:) & * WORK4(:,:,1) ) SF_SLY(:,:,2,kk,k) = reference_depth(kk) * WORK5 & * ( 2.0d0 * WORK3(:,:,2) + TLT%THICKNESS(:,:) & * WORK4(:,:,2) ) endwhere !----------------------------------------------------------------------- ! ! transition layer: use quadratic interpolation (in streamfunction) ! !----------------------------------------------------------------------- where ( reference_depth(kk) > TLT%DIABATIC_DEPTH(:,:) & .and. reference_depth(kk) <= TLT%INTERIOR_DEPTH(:,:) & .and. k <= KMT(:,:) ) WORK7 = (TLT%DIABATIC_DEPTH(:,:) & - reference_depth(kk))**2 SF_SLX(:,:,1,kk,k) = - WORK7 * WORK6 & * ( WORK1(:,:,1) + TLT%INTERIOR_DEPTH(:,:) & * WORK2(:,:,1) ) & + reference_depth(kk) * WORK5 & * ( 2.0d0 * WORK1(:,:,1) + TLT%THICKNESS(:,:) & * WORK2(:,:,1) ) SF_SLX(:,:,2,kk,k) = - WORK7 * WORK6 & * ( WORK1(:,:,2) + TLT%INTERIOR_DEPTH(:,:) & * WORK2(:,:,2) ) & + reference_depth(kk) * WORK5 & * ( 2.0d0 * WORK1(:,:,2) + TLT%THICKNESS(:,:) & * WORK2(:,:,2) ) SF_SLY(:,:,1,kk,k) = - WORK7 * WORK6 & * ( WORK3(:,:,1) + TLT%INTERIOR_DEPTH(:,:) & * WORK4(:,:,1) ) & + reference_depth(kk) * WORK5 & * ( 2.0d0 * WORK3(:,:,1) + TLT%THICKNESS(:,:) & * WORK4(:,:,1) ) SF_SLY(:,:,2,kk,k) = - WORK7 * WORK6 & * ( WORK3(:,:,2) + TLT%INTERIOR_DEPTH(:,:) & * WORK4(:,:,2) ) & + reference_depth(kk) * WORK5 & * ( 2.0d0 * WORK3(:,:,2) + TLT%THICKNESS(:,:) & * WORK4(:,:,2) ) endwhere !----------------------------------------------------------------------- ! ! interior, adiabatic region: no interpolation is needed. note that ! "dzw" is introduced here, too, for consistency. ! !----------------------------------------------------------------------- where ( reference_depth(kk) > TLT%INTERIOR_DEPTH(:,:) & .and. k <= KMT(:,:) ) SF_SLX(:,:,1,kk,k) = KAPPA_THIC(:,:,kk,k) & * SLX(:,:,1,kk,k) * (1.0d0/odz(k)) SF_SLX(:,:,2,kk,k) = KAPPA_THIC(:,:,kk,k) & * SLX(:,:,2,kk,k) * (1.0d0/odz(k)) SF_SLY(:,:,1,kk,k) = KAPPA_THIC(:,:,kk,k) & * SLY(:,:,1,kk,k) * (1.0d0/odz(k)) SF_SLY(:,:,2,kk,k) = KAPPA_THIC(:,:,kk,k) & * SLY(:,:,2,kk,k) * (1.0d0/odz(k)) endwhere enddo ! end of kk-loop enddo ! end of k-loop !----------------------------------------------------------------------- !EOC end subroutine merged_streamfunction !*********************************************************************** !BOP ! !IROUTINE: apply_vertical_profile_to_isop_hor_diff ! !INTERFACE: subroutine apply_vertical_profile_to_isop_hor_diff ! !DESCRIPTION: ! Apply vertical tapers to KAPPA_ISOP and HOR_DIFF based on their ! vertical location with respect to the diabatic, transition, and ! adiabatic regions. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer :: k, kk ! loop indices real*8, dimension(2) :: reference_depth !----------------------------------------------------------------------- ! ! start of tapering ! !----------------------------------------------------------------------- do k=1,k1 reference_depth(ktp) = (z(2*k) - z(1)) - 0.25d0 * (1.0d0/odz(k)) reference_depth(kbt) = (z(2*k) - z(1)) + 0.25d0 * (1.0d0/odz(k)) do kk=ktp,kbt !----------------------------------------------------------------------- ! ! diabatic region: no isopycnal diffusion ! !----------------------------------------------------------------------- where ( reference_depth(kk) <= TLT%DIABATIC_DEPTH(:,:) & .and. k <= KMT(:,:) ) KAPPA_ISOP(:,:,kk,k) = 0.d0 endwhere !----------------------------------------------------------------------- ! ! transition layer: a linear combination of isopcynal and horizontal ! diffusion coefficients ! !----------------------------------------------------------------------- where ( reference_depth(kk) > TLT%DIABATIC_DEPTH(:,:) & .and. reference_depth(kk) <= TLT%INTERIOR_DEPTH(:,:) & .and. k <= KMT(:,:) .and. & TLT%THICKNESS(:,:) > eps ) HOR_DIFF(:,:,kk,k) = ( TLT%INTERIOR_DEPTH(:,:) & - reference_depth(kk) ) * HOR_DIFF(:,:,kk,k) & / TLT%THICKNESS(:,:) KAPPA_ISOP(:,:,kk,k) = ( reference_depth(kk) & - TLT%DIABATIC_DEPTH(:,:) ) & * KAPPA_ISOP(:,:,kk,k) / TLT%THICKNESS(:,:) endwhere !----------------------------------------------------------------------- ! ! interior region: no horizontal diffusion ! !----------------------------------------------------------------------- where ( reference_depth(kk) > TLT%INTERIOR_DEPTH(:,:) & .and. k <= KMT(:,:) ) HOR_DIFF(:,:,kk,k) = 0.d0 endwhere enddo ! end of kk-loop enddo ! end of k-loop !----------------------------------------------------------------------- !EOC end subroutine apply_vertical_profile_to_isop_hor_diff !*********************************************************************** end module hmix_GM