!/ Origination: Apr 2015 !/ Author : Brandon G. Reichl !/ MODULE KPP USE setvars IMPLICIT NONE ! These must be public for the model to work public :: LASTX, LASTY, KPPH, KPP_Lagrangian_LOG ! These must be public for expaned KPP outputs public :: LA, KMKPP, KHKPP, KPPRIB, NN, HMO, KPPb private !Switches for kpp.nml ! KPP_LT: Switch for enhancement due to Langmuir turbulence LOGICAL :: KPP_LT_LOG ! KPP_LAGRANGIAN: Switch for using Lagrangian currents in KPP LOGICAL :: KPP_LAGRANGIAN_LOG ! KPP_SIG_MATCHING: Switch for matching mixing at base to interior LOGICAL :: KPP_SIG_MATCH_SBL_LOG LOGICAL :: KPP_SIG_MATCH_BBL_LOG LOGICAL :: KPP_SBL_LOG LOGICAL :: KPP_BBL_LOG LOGICAL :: KPP_INT_LOG ! KPP_STD_Out: Switch for diagnostic KPP output to STD out. LOGICAL :: KPP_STD_Out_LOG ! KPP_HEATLFUX_LOG: Switch to disable convective mixing within KPP SBL LOGICAL :: KPP_HEATFLUX_LOG ! USE_KPP_K: Switch to put KPP into diagnostic mode LOGICAL :: USE_KPP_K_LOG ! Limit KPP based on Ekman depth LOGICAL :: KPP_SBL_Ek_LIM_LOG REAL :: MIN_KPP_KH, MIN_KPP_KM REAL :: MAX_KPP_KH, MAX_KPP_KM REAL :: KPP_RIc ! Shared parameters for KPP_ ! Moving here from pom.h to make easier for switching b/t coupling modes REAL, DIMENSION(:,:), ALLOCATABLE :: KPPb, KPPh, LA, LASTX, $ LASTY, USTAR2, HEK, HLim, $ HMO, BR, $ BFLUX, HFlux, SFlux REAL, DIMENSION(:,:,:), ALLOCATABLE :: khkpp, kmkpp, kpprib, $ SS, NN INTEGER, DIMENSION(:,:), ALLOCATABLE :: KPPhind ! Similarity coefficients real, parameter:: zetas = -1.0 real, parameter:: zetam = -0.2 real, parameter:: cs = 98.96 real, parameter:: Cm = 8.38 real, parameter:: As = -28.86 real, parameter:: Am = 1.26 ! For non-local flux (Cs differs from cs, naming 'See-es') ! Solved for from Cs= kappa * cstar * (kappa * epsilon * cs)^(1/3) real, parameter:: SEEes = 6.3275 ! Heat capacity of water [J/kg/degC] real, parameter:: Cp = 3985 ! Critical gradient Richardson number for turbulent mixing real, parameter :: Ri0 = 0.7 ! Max interior viscosity/diffusivity due to shear instability real, parameter :: KMRIg = 10.0E-4 real, parameter :: KHRIg = 10.0E-4 ! Interior viscosity due to internal wave breaking real, parameter :: KMWB = 1.e-5 real, parameter :: KHWB = 1.e-6 ! buoyancy frequency (1/s2) limit for convection real, parameter :: bvfcon = -2.e-5 ! Interior viscosity due to convective diffusivity due to shear instability real, parameter :: KC = 0.01 ! Hard-coded KPP coefficients (note changing these may require ! modification in WW3 as well) REAL, PARAMETER :: KPP_SL_DPTnd = 0.1 REAL, PARAMETER :: KPP_LASL_DPTnd = 0.2 INTEGER :: IUP, JUP ! integer, parameter :: KPPLAGRANGIAN=0 ! Should be set in input file ! public:: ProfKPP, initKPP contains Subroutine profKPP USE LINEAR_MD !Interface to KPP integer :: indi, indj, indk REAL :: ALPHA_TE !Thermal Expansion Coefficient [1/degC] REAL :: BETA_HE REAL :: STKM real dpt(1), botu(1),botv(1),DU,DV,SHEAR_DIR,WAVE_DIR, $ SHEAR_WAVE_ANG real r(5),ad1(5),ad2(5) REAL TOUTTEST ! irradiance parameters after Paulson and Simpson (1977) ! ntp 1 2 3 4 5 ! Jerlov type i ia ib ii iii data r / .58e0, .62e0, .67e0, .77e0, .78e0 / data ad1 / .35e0, .60e0, 1.0e0, 1.5e0, 1.4e0 / data ad2 / 23.e0, 20.e0, 17.e0, 14.e0, 7.9e0 / KPPHIND=1 KMKPP=0.0 KHKPP=0.0 KPPRIB = 0.0 BFLUX = 0.0 HFLUX = 0.0 SFLUX = 0.0 if (ionedim.eq.0) then IUP = imm1 JUP = jmm1 else IUP = IM JUP = JM endif !/---------------------------------------------------------------------+ ! Get u*, B*, hek, and hmon !\---------------------------------------------------------------------+ if (ionedim.eq.0) then do indj=1,JUP do indi=1,IUP ustar2(indi,indj) = sqrt((.5e0*(wusurf(indi,indj) $ * dum(indi,indj) + wusurf(indi+1,indj) $ * dum(indi+1,indj)))**2 + $ (.5e0*(wvsurf(indi,indj)*dvm(indi,indj) $ +wvsurf(indi,indj+1)*dvm(indi,indj+1)))**2) enddo enddo call exchange2d_mpi(ustar2,im,jm) else do indj=1,JUP do indi=1,IUP ustar2(indi,indj) = sqrt((wusurf(indi,indj) $ * dum(indi,indj))**2 + $ (wvsurf(indi,indj)*dvm(indi,indj))**2) enddo enddo endif do indj=1,JUP do indi=1,IUP hek(indi,indj) = 0.7 * sqrt(ustar2(indi,indj)) $ / max(SMALL,abs(cor(indi,indj))) ! Get Total surface heat flux in W/m^2 if (KPP_HEATFLUX_LOG) then HFLUX(indi,indj)=(wtsurf(indi,indj)-swrad(indi,indj)) $ *rhoref*cp else HFLUX(indi,indj) = 0.0 endif !Get Alpha_te (thermal expansion coefficient) in 1/degC CALL GETExpansioncoeffs(S(indi,indj,1),T(indi,indj,1), $ alpha_te,beta_he) ! this is subtracted from surface flux due to POM convention. BR(indi,indj)=swrad(indi,indj) $ *(1.-( r(ntp) *exp(-kpph(indi,indj)/ad1(ntp) ) $ + (1.e0-r(ntp))*exp(-kpph(indi,indj)/ad2(ntp)))) $ *rhoref*Cp ! Get Bf in m^2/s^3 Bflux(indi,indj) = (HFLUX(indi,indj)+BR(indi,indj))* $ alpha_te * grav/rhoref/Cp ! $ + W(indi,indj,1)*Rhoref * beta_he ! $ * (S(indi,indj,1)-0.0)*grav/rhoref; ! Get Minon-Obhukov length in m if (abs(hflux(indi,indj)).ge.0.001) then hmo(indi,indj) = 1.0 * sqrt(ustar2(indi,indj))* $ ustar2(indi,indj) / kappa / BFLUX(indi,indj) else hmo(indi,indj) =sign(1.0e10,HFLUX(indi,indj)) endif enddo enddo if (ionedim.eq.0) then call exchange2d_mpi(hek,im,jm) call exchange2d_mpi(hmo,im,jm) call exchange2d_mpi(hflux,im,jm) call exchange2d_mpi(bflux,im,jm) call exchange2d_mpi(BR,im,jm) endif if (KPP_SBL_Ek_LIM_LOG) then hlim = hek else hlim = H endif !/---------------------------------------+ ! Calculate Shear and buoyancy frequency !\---------------------------------------+ call Freq_Prof() !/-----------------------------+ !Compute interior diffusivities !\-----------------------------+ IF (KPP_INT_LOG) THEN call interiorK() ENDIF !/---------------------------------------------------------------------+ ! Construct Stokes current onto POM vertical grid !\---------------------------------------------------------------------+ LA(:,:)=1.e8 US(:,:,:)=0.0 VS(:,:,:)=0.0 do indj=1,JUP do indi=1,IUP if (KPP_LT_LOG) then STKM=sqrt(LASTX(indi,indj)**2+LASTY(indi,indj)**2) if (STKM.gt.0.0001) then DPT(1) = KPPh(indi,indj)*KPP_LASL_DPTND/h(indi,indj) call linear(-999.,-999.,u(indi,indj,:), $ BOTU(1),-zz,DPT(1),kb,1) call linear(-999.,-999.,v(indi,indj,:), $ BOTV(1),-zz,DPT(1),kb,1) DU = U(indi,indj,1)-BOTU(1) DV = V(indi,indj,1)-BOTV(1) ! If not in 1d (a-grid) need to do averaging ! if (ionedim.eq.0) then call linear(-999.,-999.,u(indi+1,indj,:), $ BOTU(1),-zz,DPT(1),kb,1) call linear(-999.,-999.,v(indi,indj+1,:), $ BOTV(1),-zz,DPT(1),kb,1) DU =DU+U(indi+1,indj,1)-BOTU(1) DV =DV+V(indi,indj+1,1)-BOTV(1) DU = DU * 0.5 DV = DV * 0.5 endif SHEAR_DIR = atan2(DV,DU) WAVE_DIR = atan2(.5*(LASTY(indi,indj) $ +LASTY(indi,indj+1)), $ .5*(LASTX(indi,indj)+LASTX(indi+1,indj))) SHEAR_WAVE_ANG = max(0.0001,cos(SHEAR_DIR-WAVE_DIR)) LA(indi,indj) = max(0.1, $ sqrt(max(0.001,sqrt(ustar2(indi,indj)))/STKM)) LA(indi,indj) = LA(indi,indj)*min(10.0, $ (sqrt(1./SHEAR_WAVE_ANG))) else LA(indi,indj) = 1.e10 endif else LA(indi,indj)=1.e10 endif enddo enddo if (ionedim.eq.0) then call exchange2d_mpi(LA,im,jm) endif !/---------------------------------------------------------------------+ ! Interface to SBL routines !\---------------------------------------------------------------------+ if (KPP_SBL_LOG) then call ComputeRIb() call GetKPPH() call KPPsurface() else KPPRIB=0.0 KPPH = 0.0 endif !/---------------------------------------------------------------------+ ! Interface to BBL routines !\---------------------------------------------------------------------+ if (KPP_BBL_LOG) then call KPPbottom endif !/---------------------------------------------------------------------+ ! Apply KPP !\---------------------------------------------------------------------+ if (use_KPP_K_LOG) THEN DO INDI=1,IM DO INDJ=1,JM DO INDK=1,KB km(indi,indj,indk)=min(max(Min_KPP_KM, $ kmkpp(indi,indj,indk)),Max_KPP_KM) kh(indi,indj,indk)=min(max(Min_KPP_KH, $ khkpp(indi,indj,indk)),Max_KPP_KH) ENDDO ENDDO ENDDO endif !/---------------------------------------------------------------------+ ! Output Diagnostics !\---------------------------------------------------------------------+ if (KPP_STD_OUT_LOG) THEN write(*,*)'UST2:',minval(USTAR2),maxval(USTAR2) write(*,*)'HEAT:',minval(HFLUX),maxval(HFLUX) write(*,*)'Buoy:',minval(Bflux),maxval(BFLUX) write(*,*)'STK: ',maxval(LASTX),maxval(LASTY) write(*,*)'LA: ',minval(LA),maxval(LA) write(*,*)'SBL: ',minval(KPPH),maxval(KPPH) write(*,*)'KM: ',minval(KMKPP),maxval(KMKPP) write(*,*)'KH: ',minval(KHKPP),maxval(KHKPP) endif ! cosmetics: make boundr. values as interior (even if not used, printout ! may show strange values) do indk=1,kb do indi=1,im if(n_north.eq.-1) then km(indi,jm,indk)=km(indi,jmm1,indk) kh(indi,jm,indk)=kh(indi,jmm1,indk) end if if(n_south.eq.-1) then km(indi,1,indk)=km(indi,2,indk) kh(indi,1,indk)=kh(indi,2,indk) end if end do do indj=1,jm if(n_east.eq.-1) then km(im,indj,indk)=km(imm1,indj,indk) kh(im,indj,indk)=kh(imm1,indj,indk) end if if(n_west.eq.-1) then km(1,indj,indk)=km(2,indj,indk) kh(1,indj,indk)=kh(2,indj,indk) end if end do end do do indk=1,kb do indi=1,im do indj=1,jm km(indi,indj,indk)=km(indi,indj,indk)*fsm(indi,indj) kh(indi,indj,indk)=kh(indi,indj,indk)*fsm(indi,indj) end do end do end do RETURN end Subroutine profKPP !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- SUBROUTINE initKPP namelist/KPP_nml/ $ KPP_LT_LOG, $ KPP_LAGRANGIAN_LOG, $ KPP_SIG_MATCH_SBL_LOG, $ KPP_SIG_MATCH_BBL_LOG, $ KPP_SBL_LOG, $ KPP_INT_LOG, $ KPP_BBL_LOG, $ KPP_STD_Out_LOG, $ KPP_HEATFLUX_LOG, $ USE_KPP_K_LOG, $ KPP_SBL_Ek_LIM_LOG, $ MIN_KPP_KH, $ MIN_KPP_KM, $ MAX_KPP_KH, $ MAX_KPP_KM, $ KPP_RIc ! read input namelist open(74,file='kpp.nml',status='old') read(74,nml=KPP_nml) close(74) if (KPP_LAGRANGIAN_LOG) then print*,'FATAL: Lagrangian mixing not coded in 3-way' print*,' system because Stokes drift not passed.' print*,'****************' print*,'* Quitting now *' print*,'****************' stop endif !Allocate 2-d arrays ALLOCATE(KPPb(IM,JM)) ; KPPB = 0.0; ALLOCATE(KPPH(IM,JM)) ; KPPH = 0.0; ALLOCATE(LA(IM,JM)) ; LA = 0.0; ALLOCATE(LASTX(IM,JM)) ; LASTX = 0.0; ALLOCATE(LASTY(IM,JM)) ; LASTY = 0.0; ALLOCATE(USTAR2(IM,JM)); USTAR2 = 0.0; ALLOCATE(HEK(IM,JM)) ; HEK = 0.0; ALLOCATE(HLIM(IM,JM)) ; HLIM = 0.0; ALLOCATE(HMO(IM,JM)) ; HMO = 0.0; ALLOCATE(BR(IM,JM)) ; BR = 0.0; ALLOCATE(BFLUX(IM,JM)) ; BFLUX = 0.0; ALLOCATE(HFLUX(IM,JM)) ; HFLUX = 0.0; ALLOCATE(SFLUX(IM,JM)) ; SFLUX = 0.0; ALLOCATE(KPPHIND(IM,JM)); KPPHIND = 0; !Allocate 3-d arrays ALLOCATE(KMKPP(IM,JM,KB)) ; KMKPP = 0.0; ALLOCATE(KHKPP(IM,JM,KB)) ; KHKPP = 0.0; ALLOCATE(KPPRIB(IM,JM,KB)); KPPRIB = 0.0; ALLOCATE(SS(IM,JM,KB)) ; SS = 0.0; ALLOCATE(NN(IM,JM,KB)) ; NN = 0.0; return END SUBROUTINE initKPP !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- SUBROUTINE getexpansioncoeffs(S,T,ALPHA_TE,BETA_HE) !--------------------------------------------------- ! This subroutine computes the thermal expansion ! coefficient at a given S/T ! Coded into POM following GOTM by Brandon G. Reichl ! Inputs: S, T ! S: Salinity [ppt] ! T: Temperature [deg C] ! Outputs: ALPHA_TE, BETA_HE ! ALPHA_TE: Thermal expansion [1/deg C] ! BETA_TE: Haline expansion [1/ppt] !--------------------------------------------------- real, intent(in) :: S, T real, intent(out) :: ALPHA_TE,BETA_HE real :: Buoy1, Buoy2 real*8 :: Rho1, Rho2 REAL, parameter :: DELTA=0.01 !------------------------------------- !1. Get Alpha call DPeqstate(S,T+0.5*DELTA,0.0,Rho1) call DPeqstate(S,T-0.5*DELTA,0.0,Rho2) Buoy1 = (Rho1)/rhoref Buoy2 = (Rho2)/rhoref ALPHA_TE = (Buoy1 - Buoy2) / DELTA !------------------------------------- !2. Get Beta call DPeqstate(S+0.5*DELTA,T,0.0,Rho1) call DPeqstate(S-0.5*DELTA,T,0.0,Rho2) Buoy1 = (Rho1)/rhoref Buoy2 = (Rho2)/rhoref BETA_HE = (Buoy1 - Buoy2) / DELTA !------------------------------------- return end Subroutine getexpansioncoeffs !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- SUBROUTINE eqstate(SI,TI,DPT, RHO) !--------------------------------------------------------- ! This subroutines does the identical equation of state ! to existing POM routine. Copied here since the version ! is not a module, if state equation is changed this ! must be considered. ! Inputs: SI, TI, DPT ! SI: Salinity ! TI: Temperature ! DPT: Depth (used for pressure, set to 0 for ! computing potential density) ! Output: RHO: Density solving Unesco eq. of state !-- ! Shared variables accessed ! grav - pom.h ! rhoref - pom.h !--------------------------------------------------------- !/ REAL, INTENT(IN) :: TI, SI,DPT REAL, INTENT(OUT) :: RHO !/ REAL :: cr, tr,sr,tr2,tr3,tr4,p,rhor !/ tr=ti sr=si tr2=tr*tr tr3=tr2*tr tr4=tr3*tr p=grav*rhoref*(DPT)*1.e-5 rhor=-0.157406e0 $ +6.793952e-2*tr $ -9.095290e-3*tr2 $ +1.001685e-4*tr3 $ -1.120083e-6*tr4 $ +6.536332e-9*tr4*tr rhor=rhor+(0.824493e0-4.0899e-3*tr $ +7.6438e-5*tr2-8.2467e-7*tr3 $ +5.3875e-9*tr4)*sr $ +(-5.72466e-3+1.0227e-4*tr $ -1.6546e-6*tr2)*abs(sr)**1.5 $ +4.8314e-4*sr*sr cr=1449.1e0+.0821e0*p+4.55e0*tr-.045e0*tr2 $ +1.34e0*(sr-35.e0) rhor=rhor+1.e5*p/(cr*cr)*(1.e0-2.e0*p/(cr*cr)) RHO=rhor+1000.e0 RETURN end SUBROUTINE eqstate !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- SUBROUTINE DPeqstate(SIr,TIr,DPTr,rho) !--------------------------------------------------------- ! ***SIMILAR TO ABOVE BUT COMPUTES AT DOUBLE PRECISION*** ! ***AND RETURNS THE DENSITY MINUS 1,000*** !--------------------------------------------------------- ! This subroutines does the identical equation of state ! to existing POM routine. Copied here since the version ! is not a module, if state equation is changed this ! must be considered. ! Inputs: SI, TI, DPT ! SI: Salinity ! TI: Temperature ! DPT: Depth (used for pressure, set to 0 for ! computing potential density) ! Output: RHO: Density solving Unesco eq. of state !-- ! Shared variables accessed ! grav - pom.h ! rhoref - pom.h !--------------------------------------------------------- !/ real, intent (IN) :: TIr, SIr,DPTr real*8, INTENT(OUT) :: rho real*8 cr, tr,sr,tr2,tr3,tr4,p,rhor real*8 TI, SI, DEPTH TI = TIr*1.d+0 SI = SIr*1.d+0 DEPTH = DPTr*1.d+0 tr=ti sr=si tr2=tr*tr tr3=tr2*tr tr4=tr3*tr p=grav*rhoref*(DEPTH)*1.d-5 rhor=-0.157406d0+6.793952d-2*tr $ -9.095290d-3*tr2+1.001685d-4*tr3 $ -1.120083d-6*tr4+6.536332d-9*tr4*tr rhor=rhor+(0.824493d0-4.0899d-3*tr $ +7.6438d-5*tr2-8.2467d-7*tr3 $ +5.3875d-9*tr4)*sr $ +(-5.72466d-3+1.0227d-4*tr $ -1.6546d-6*tr2)*abs(sr)**1.5 $ +4.8314d-4*sr*sr cr=1449.1d0+.0821d0*p+4.55d0*tr-.045d0*tr2 $ +1.34d0*(sr-35.d0) rhor=rhor+1.d5*p/(cr*cr)*(1.d0-2.d0*p/(cr*cr)) RHO=rhor RETURN end SUBROUTINE DPeqstate !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- SUBROUTINE VSCALE(BF,UST,DPT,WM,WS) !-----------------------------------------------------+ ! This subroutine calculates the velocity scale ! used within KPP based on the surface Buoyancy ! flux, the surface momentum flux, and the depth ! Input: BF, UST, DPT ! BF: Surface Buoyancy flux ! UST: Surface friction velocity ! DPT: Depth ! OUTPUT: WM, WS ! WM: Velocity scale for momentum ! WS: Velocity scale for scalars ! Adapted into POM following GOTM by Brandon G. Reichl !-----------------------------------------------------+ REAL, intent(in):: UST, DPT, BF REAL, intent(out) :: WM, WS REAL :: ZETA, ZETAH ZETAH = kappa * DPT * BF ZETA = ZETAH / (UST**3+1.e-20); if (ZETAH.ge.0.) then !/Stable WM=kappa*UST/(1.+5.*ZETA); WS=WM; else !/Unstable if (ZETA.gt.ZETAm) then WM=kappa*UST*((1.-16.*ZETA)**(.25)) else WM=kappa*((UST**3*Am-Cm*ZETAH)**(1./3.)) endif if (ZETA.gt.ZETAs) then WS=kappa*UST*((1.-16.*ZETA)**(.5)) else WS=kappa*((UST**3*As-Cs*ZETAH)**(1./3.)) endif endif return end subroutine Vscale !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- Subroutine Freq_Prof() !----------------------------------------------------------- ! This subroutine computes the profiles of the buoyancy and ! shear frequencies. ! Brandon G. Reichl !-- ! Shared variables accessed ! SS - KPP ! NN - KPP ! grav - pom.h ! rhoref - pom.h ! ionedim - pom.h ! fsm(im,jm) - pom.h ! T(im,jm,kb) - pom.h ! S(im,jm,kb) - pom.h ! dzz(kb) - pom.h ! h(im,jm) - pom.h !----------------------------------------------------------- integer :: indj, indi, indk real :: depth, uup, vup, udn, vdn real*8 :: rhou, rhol, drhoz !-------------------------------- SS = 0.0 NN = 1.d-8 do indj=1,JUP do indi=1,IUP if (fsm(indi,indj).gt.0.1) then SS(indi,indj,:)=0. NN(indi,indj,:)=0.0!1.e-8 do indk=2,kbm1 ! Depth of Mid-Point (computes NN at depth) DEPTH = 0.0 - zz(indk) * h(indi,indj) call DPeqstate(S(indi,indj,indk-1), $ T(indi,indj,indk),DEPTH,RHOU) call DPeqstate(S(indi,indj,indk), $ T(indi,indj,indk), $ DEPTH,RHOL) DRHOz = (RHOU - RHOL) call DPeqstate(S(indi,indj,indk), $ T(indi,indj,indk-1),DEPTH,RHOU) call DPeqstate(S(indi,indj,indk), $ T(indi,indj,indk), $ DEPTH,RHOL) DRHOz = DRHOZ + (RHOU - RHOL) NN(indi,indj,indk) = ((-grav/RHOREF $ *(DRHOz)/((dzz(indk-1))*h(indi,indj)))) if (ionedim.eq.0) then Uup=.5*(u(indi,indj,indk-1)*dum(indi,indj) $ +u(indi+1,indj,indk-1)*dum(indi+1,indj)) Vup=.5*(v(indi,indj,indk-1)*dvm(indi,indj) $ +v(indi,indj+1,indk-1)*dvm(indi,indj+1)) Udn=.5*(u(indi,indj,indk)*dum(indi,indj) $ +u(indi+1,indj,indk)*dum(indi+1,indj)) Vdn=.5*(v(indi,indj,indk)*dvm(indi,indj) $ +v(indi,indj+1,indk)*dvm(indi,indj+1)) else Uup = u(indi,indj,indk-1)*dum(indi,indj) Udn = u(indi,indj,indk)*dum(indi,indj) Vup = v(indi,indj,indk-1)*dvm(indi,indj) Vdn = v(indi,indj,indk)*dvm(indi,indj) endif SS(indi,indj,indk) = ((Uup-Udn)**2 +(Vup-Vdn)**2) $ /(dzz(indk)*h(indi,indj))**2 enddo endif enddo enddo if (ionedim.eq.0) then do indk=1,kb call exchange2d_mpi(SS(:,:,indk),im,jm) call exchange2d_mpi(NN(:,:,indk),im,jm) enddo endif return end subroutine Freq_Prof !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- Subroutine InteriorK() !-------------------------------------------------------- ! This subroutine computes the interior viscosity for the ! KPP scheme. ! Adapted into POM following GOTM by Brandon G. Reichl !-- ! Shared variables accessed ! jmm1, imm1 - pom.h ! fs (indi,indj) - pom.h ! NN - KPP ! SS - KPP ! Ri0 - KPP ! KC - KPP ! KMWB - KPP ! KHWB - KPP ! KMRIg - KPP ! KHRIg - KPP ! bvcon - KPP !--------------------------------------------------------- real :: RIG(kb) real*8 :: nu_sx, nu_sxc, shear2 integer :: indj, indi, indk,iol do indj=1,JUP do indi=1,IUP if (fsm(indi,indj).gt.0.1) then do indk=1,kb RIG(indk)=NN(indi,indj,indk) $ /(SS(indi,indj,indk)+1.E-14) enddo do indk=kbm1,1,-1 if (indk.ge.2.and.indk.lt.kb) then ! Gradient Richardson mixing (smoothing) RIG(indk)=0.25*Rig(indk-1) + 0.50*Rig(indk) $ + 0.25*Rig(indk+1) nu_sx = min(1.,max(0.,Rig(indk))/Ri0)*1.d+0 nu_sx = (1.0 - nu_sx*nu_sx) nu_sx = nu_sx * nu_sx * nu_sx shear2 = SS(indi,indj,indk)*1.d+0 nu_sx=nu_sx*shear2*shear2/ $ (shear2*shear2+16.0E-10) endif nu_sxc = max(NN(indi,indj,indk),bvfcon)*1.d+0 nu_sxc = min(1.,(bvfcon-nu_sxc)/bvfcon) nu_sxc = 1.d+0 - nu_sxc*nu_sxc nu_sxc = nu_sxc*nu_sxc*nu_sxc kmkpp(indi,indj,indk)= real(nu_sxc*KC) $ +real(nu_sx*KMRIg) + KMWB khkpp(indi,indj,indk)= real(nu_sxc*KC) $ +real(nu_sx*KHRIg) + KHWB enddo endif enddo enddo if (ionedim.eq.0) then do indk=1,kb call exchange2d_mpi(KMKPP(:,:,indk),im,jm) call exchange2d_mpi(KHKPP(:,:,indk),im,jm) enddo endif return end subroutine InteriorK !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- Subroutine ComputeRIb() !-------------------------------------- ! This subroutine computes the Surface layer bulk Richardson ! number for use in the KPP SBL model. ! Brandon G. Reichl !-- ! Share variables accessed ! jmm1,imm1 - pom.h ! fsm(im,jm) - pom.h ! dum(im,jm) - pom.h ! dbm(im,jm) - pom.h ! dz(kb) - pom.h ! h(im,jm) - pom.h ! S(im,jm,kb) -pom.h ! T(im,jm,kb) - pom.h ! U(im,jm,kb) - pom.h ! US(im,jm,kb) - pom.h ! VS(im,jm,kb) - pom.h ! KPP_RIc - KPP ! KPPRIb(im,jm,kb) - KPP ! --------------------------------------------------- real, dimension(im,jm) :: AVGDPTH integer :: indi, indj, indk logical :: first integer :: kstrt real :: ustop, vstop, DEPTH, USURF, VSURF, DUz, DVz, $ RItop, FLT, VT2, VTC, RIBOT, WM, WS real*8 :: RHOU,RHOL,DRHOZ,RHOLEV,RHOSURF !/ VTC = 1.6*sqrt(.2)/(sqrt(98.96*0.1)*KPP_Ric*.4*.4) !/ do indj=1,JUP do indi=1,IUP if (fsm(indi,indj).gt.0.1) then first=.true. kstrt=1 ustop=0.0; vstop=0.0 ! Determine thickness of top level DEPTH = dz(1) * h(indi,indj) ! compute density of top level call DPeqstate(S(indi,indj,1),T(indi,indj,1), $ 0.0,rhosurf) rhosurf=(rhosurf)*DEPTH if (kpp_lagrangian_log) then if (ionedim.eq.0) then ustop = 0.5e0*(us(indi,indj,1)*dum(indi,indj) $ +us(indi+1,indj,1)*dum(indi+1,indj)) vstop = 0.5e0*(vs(indi,indj,1)*dvm(indi,indj) $ +vs(indi,indj+1,1)*dvm(indi,indj+1)) else ustop = (us(indi,indj,1)*dum(indi,indj)) vstop = (vs(indi,indj,1)*dvm(indi,indj)) endif endif if (ionedim.eq.0) then usurf = (.5e0*(u(indi,indj,1)*dum(indi,indj) $ +u(indi+1,indj,1)*dum(indi+1,indj))+ustop) $ *DEPTH vsurf = (.5e0*(v(indi,indj,1)*dvm(indi,indj) $ +v(indi,indj+1,1)*dvm(indi,indj+1))+vstop) $ *DEPTH else usurf = (u(indi,indj,1)*dum(indi,indj) + ustop)*DEPTH vsurf = (v(indi,indj,1)*dvm(indi,indj) + vstop)*DEPTH endif ! Determine how deep to average for surface layer referenced value AVGDPTH(indi,indj) = $ max(KPPh(indi,indj)*KPP_SL_DPTnd,DEPTH) do indk=2,kbm1 !Depth of bottom of cell to see if whole cell is ! within SBL DEPTH = 0.0 - z(indk+1) * h(indi,indj) if (depth.lt.AVGDPTH(indi,indj)) then ! Getting u-component/v-component if (kpp_lagrangian_log) then if (ionedim.eq.0) then ustop = 0.5e0*( $ us(indi,indj,indk)*dum(indi,indj) $ +us(indi+1,indj,indk)*dum(indi+1,indj)) vstop = 0.5e0*( $ vs(indi,indj,indk)*dvm(indi,indj) $ +vs(indi,indj+1,indk)*dvm(indi,indj+1)) else ustop = us(indi,indj,indk)*dum(indi,indj) vstop = vs(indi,indj,indk)*dvm(indi,indj) endif endif if (ionedim.eq.0) then usurf = usurf+(.5e0* $ (u(indi,indj,indk)*dum(indi,indj) $ +u(indi+1,indj,indk)*dum(indi+1,indj)) $ +ustop)*dz(indk)*h(indi,indj) vsurf = vsurf+(.5e0* $ (v(indi,indj,indk)*dvm(indi,indj) $ +v(indi,indj+1,indk)*dvm(indi,indj+1)) $ +vstop)*dz(indk)*h(indi,indj) else usurf = usurf+(u(indi,indj,indk)*dum(indi,indj) $ +ustop)*dz(indk)*h(indi,indj) vsurf = vsurf+(v(indi,indj,indk)*dvm(indi,indj) $ +vstop)*dz(indk)*h(indi,indj) endif call DPeqstate(S(indi,indj,indk),T(indi,indj,indk), $ 0.0,rholev) rhosurf = rhosurf + $ (rholev)*dz(indk)*h(indi,indj) elseif (first) then kstrt=indk if (kpp_lagrangian_log) then if (ionedim.eq.0) then ustop = 0.5e0*( $ us(indi,indj,indk)*dum(indi,indj) $ +us(indi+1,indj,indk)*dum(indi+1,indj)) vstop = 0.5e0*( $ vs(indi,indj,indk)*dvm(indi,indj) $ +vs(indi,indj+1,indk)*dvm(indi,indj+1)) else ustop = us(indi,indj,indk)*dum(indi,indj) vstop = vs(indi,indj,indk)*dvm(indi,indj) endif endif if (ionedim.eq.0) then usurf = usurf+(.5e0* $ (u(indi,indj,indk)*dum(indi,indj) $ +u(indi+1,indj,indk)*dum(indi+1,indj)) $ +ustop)*(AVGDPTH(indi,indj) $ +z(indk)*h(indi,indj)) vsurf = vsurf+(.5e0* $ (v(indi,indj,indk)*dvm(indi,indj) $ +v(indi,indj+1,indk)*dvm(indi,indj+1)) $ +vstop)*(AVGDPTH(indi,indj) $ +z(indk)*h(indi,indj)) else usurf = usurf + $ (u(indi,indj,indk)*dum(indi,indj) $ +ustop)*(AVGDPTH(indi,indj) $ +z(indk)*h(indi,indj)) vsurf = vsurf+ $ (v(indi,indj,indk)*dvm(indi,indj) $ +vstop)*(AVGDPTH(indi,indj) $ +z(indk)*h(indi,indj)) endif call DPeqstate(S(indi,indj,indk),T(indi,indj,indk), $ 0.0,rholev) rhosurf =rhosurf+(rholev)*(AVGDPTH(indi,indj) $ +z(indk)*h(indi,indj)) first=.false. endif enddo !/ usurf = usurf/AVGDPTH(indi,indj) vsurf = vsurf/AVGDPTH(indi,indj) rhosurf = rhosurf/AVGDPTH(indi,indj) !/ KPPRIB(indi,indj,:)=0.0 do indk=kstrt,kbm1 if (kpp_lagrangian_log) then if (ionedim.eq.0) then ustop = 0.5e0*( $ us(indi,indj,indk)*dum(indi,indj) $ +us(indi+1,indj,indk)*dum(indi+1,indj)) vstop = 0.5e0*( $ vs(indi,indj,indk)*dvm(indi,indj) $ +vs(indi,indj+1,indk)*dvm(indi,indj+1)) else ustop = us(indi,indj,indk)*dum(indi,indj) vstop = vs(indi,indj,indk)*dvm(indi,indj) endif endif !Depth on ZZ since u DEPTH = 0.0 - zz(indk) * h(indi,indj) !is surf 0.0? et/el? if (ionedim.eq.0) then DUz = usurf - (.5e0*( $ u(indi,indj,indk)*dum(indi,indj) $ +u(indi+1,indj,indk)*dum(indi+1,indj))+ $ ustop) DVz = vsurf - (.5e0*( $ v(indi,indj,indk)*dvm(indi,indj) $ +v(indi,indj+1,indk)*dvm(indi,indj+1))+ $ vstop) else DUz = usurf - ( u(indi,indj,indk)*dum(indi,indj) $ + ustop ) DVz = vsurf - ( v(indi,indj,indk)*dvm(indi,indj) $ + vstop ) endif call DPeqstate(S(indi,indj,indk), $ T(indi,indj,indk),0.0,RHOU) DRHOz=min(0.0,rhosurf-(RHOU)) RItop = -grav / rhoref * DRHOz * DEPTH FLT = min(1.0+2.3/(LA(indi,indj)**.5),50.0) if (Bflux(indi,indj).ge.0.0) then call VSCALE(Bflux(indi,indj), $ sqrt(ustar2(indi,indj)),DEPTH,WM,WS) else call VSCALE(Bflux(indi,indj), $ sqrt(ustar2(indi,indj)), $ min(0.1*KPPH(indi,indj),DEPTH),WM,WS) endif VT2 = VTC * DEPTH * sqrt(abs(NN(indi,indj,indk))) $ * WS*FLT RIbot = DUz**2 + DVz**2 + VT2 KPPRIB (indi,indj,indk) = RItop / (RIbot + $ 1.e-20) enddo endif enddo enddo if (ionedim.eq.0) then do indk=1,kb call exchange2d_mpi(KPPRIB(:,:,indk),im,jm) enddo endif return end subroutine ComputeRIb !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- subroutine GetKPPh() !--------------------------------------------------------------- ! This subroutine computes the KPP surface boundary layer depth ! from the bulk Richardson number. ! Brandon G. Reichl !-- ! Shared variables accessed ! fsm - pom.h ! ZZ, DZZ, h, imm1,jmm1,kbm1 ! KPP_RIc - KPP ! KPPh(im,jm) - KPP ! KPPRIb(im,jm,kb) - KPP !--------------------------------------------------------------- integer :: indj, indi, indk logical :: TRIG kpphind=1 kpph=1.0!Okay to reset kpph here, after used in bulk RIb do indj=1,JUP do indi=1,IUP if (fsm(indi,indj).gt.0.1) then TRIG=.true. KPPh(indi,indj) = 1.0 do indk=2,kbm1 if (KPPRIB(indi,indj,indk).gt.KPP_RIc .AND.trig) then KPPh(indi,indj) = 0.0 - zz(indk) * h(indi,indj) KPPh(indi,indj) = KPPh(indi,indj)- $ (KPPRIB(indi,indj,indk)-KPP_Ric)*dzz(indk) $ *h(indi,indj)/(KPPRIB(indi,indj,indk)- $ KPPRIB(indi,indj,indk-1)) KPPH(indi,indj)=min(hlim(indi,indj), $ kpph(indi,indj)) trig=.false. endif enddo if (trig) then KPPh(indi,indj)=abs(Z(kb)*h(indi,indj)) endif ! Limit Boundary layer depth to 1 grid cell minimum if (kpph(indi,indj).le.abs(h(indi,indj)*(Z(2)))) then KPPH(indi,indj)=abs(h(indi,indj)*(Z(2))) endif KPPHind(indi,indj)=kb trig=.true. do indk=2,kb if (trig.and. $ abs(Z(indk)*h(indi,indj)).gt.KPPh(indi,indj)) then KPPHind(indi,indj)=indk trig=.false. endif enddo else KPPh(indi,indj) = 0.0 KPPhind(indi,indj) = 1 endif enddo enddo if (ionedim.eq.0) then call exchange2d_mpi(KPPhind,im,jm) call exchange2d_mpi(KPPh,im,jm) endif return end subroutine getKPPH !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- subroutine KPPsurface !------------------------------------------------- ! This subroutine applies KPP surface boundary layer ! to get mixing coefficients and non-local mixing ! profiles. ! Some parts of this code were adapted from GOTM. ! Brandon G. Reichl !---------------------------------------------------- REAL :: WM, WS, DEPTH, SIG, FLT, FLT2 INTEGER :: KSBL, indk, indi, indj REAL :: CFF, CFF_DN, CFF_UP, K_BL, dk_bl, gm1, dgm1ds, gdm1ds, $ gt1, dgt1ds, f1, a1, a2, a3, gm, gt GAMMAT=0.0 GAMMAS=0.0 GAMMAM=0.0 do indi=1,IUP do indj=1,JUP if (fsm(indi,indj).gt.0.1) then FLT = min(1.0+1./LA(indi,indj),2.25) if (KPP_SIG_MATCH_SBL_LOG) then !Compute non-dimensional shapefunction as in GOTM f1=5.*max(Bflux(indi,indj),0.)/(Ustar2(indi,indj)**2 $ +small); DEPTH=KPPh(indi,indj) if (Bflux(indi,indj).ge.0.0) then call VSCALE(Bflux(indi,indj), $ sqrt(ustar2(indi,indj)),DEPTH,WM,WS) else call VSCALE(Bflux(indi,indj), $ sqrt(ustar2(indi,indj)), $ min(0.1*KPPH(indi,indj),DEPTH),WM,WS) endif ksbl=KPPhind(indi,indj) cff=1./abs(DZ(ksbl-1)*h(indi,indj)); cff_dn=(abs(Z(ksbl)*h(indi,indj)) $ -KPPh(indi,indj))*cff cff_up=(KPPh(indi,indj)-abs(Z(ksbl-1)*h(indi,indj))) $ *cff K_bl = cff_dn*KMKPP(indi,indj,ksbl)+ $ cff_up*KMKPP(indi,indj,ksbl-1); dK_bl = cff*(KMKPP(indi,indj,ksbl) $ -KMKPP(indi,indj,ksbl-1)); Gm1 = K_bl/(abs(KPPh(indi,indj))*WM+small) dGm1dS = -dK_bl/(WM+small) + K_bl*f1; K_bl = cff_dn*KHKPP(indi,indj,ksbl) $ +cff_up*KHKPP(indi,indj,ksbl-1); dK_bl = cff*(KHKPP(indi,indj,ksbl) $ -KHKPP(indi,indj,ksbl-1)); Gt1 = K_bl/(abs(KPPh(indi,indj))*WS+small) dGt1dS = -dK_bl/(WS+small) + K_bl*f1; endif do indk=1,kb ! on z's (w levels) DEPTH = 0.0 - z(indk) * h(indi,indj) !is surf 0.0? et/el? if (DEPTH.lt.KPPH(indi,indj)) then if (Bflux(indi,indj).ge.0.0) then call VSCALE(Bflux(indi,indj), $ sqrt(ustar2(indi,indj)),DEPTH,WM,WS) else call VSCALE(Bflux(indi,indj), $ sqrt(ustar2(indi,indj)), $ min(0.1*KPPH(indi,indj),DEPTH),WM,WS) endif SIG=min(1.,DEPTH/KPPH(indi,indj)); A1 = SIG-2.0; A2 = 3.0-2.0*SIG; A3 = SIG-1.0; if (KPP_SIG_MATCH_SBL_LOG) THEN Gm = A1+A2*Gm1!+A3*dGm1dS; Gt = A1+A2*Gt1!+A3*dGt1dS; else GM=A1 GT=A1 endif FLT2=1.+(FLT-1.)*min(1.,SIG*(1.-SIG)**2.0/0.1481) kmkpp(indi,indj,indk)=DEPTH*WM*(1.+SIG*Gm)*FLT2 khkpp(indi,indj,indk)=DEPTH*WS*(1.+SIG*Gt)*FLT2 if (bflux(indi,indj).lt.0.0) then GAMMAT(indi,indj,indk)=SIG*(1-SIG)**2 $ *SEEes*max(0.0,(WTSURF(indi,indj)- $ swrad(indi,indj))+ $ BR(indi,indj)/rhoref/Cp)!subtract short wave GAMMAS(indi,indj,indk)=SIG*(1+SIG*Gm) $ *SEEes*max(0.0,WSSURF(indi,indj)) endif endif kmkpp(indi,indj,indk)=min(10.0, $ kmkpp(indi,indj,indk)) khkpp(indi,indj,indk)=min(10.0, $ khkpp(indi,indj,indk)) if (kmkpp(indi,indj,indk).lt.0.) then print*,'Less than 0 check match' print*,cff_dn/CFF,cff_up/CFF,ksbl,indk print*,Z(ksbl)*h(indi,indj),Z(ksbl-1)*h(indi,indj) $ ,KPPh(indi,indj) print*,(1+SIG*Gm),dk_bl,k_bl print*,A1,A2,GM1,A3,DGM1DS print*,cff,(KHKPP(indi,indj,ksbl), $ KHKPP(indi,indj,ksbl-1)); stop endif enddo else do indk=1,kb kmkpp(indi,indj,indk)=1.0e-5 khkpp(indi,indj,indk)=1.0e-5 enddo endif enddo enddo if (ionedim.eq.0) then do indk=1,kb call exchange2d_mpi(KHKPP(:,:,indk),im,jm) call exchange2d_mpi(KMKPP(:,:,indk),im,jm) enddo endif return end subroutine KPPsurface !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- subroutine KPPbottom ! adding bottom boundary layer for use in KPP ! note this is a simple version of the KPP boundary layer scheme ! More complexity may/may not be worth it for this application... real, dimension(im,jm) :: hek, ustarb2 real, dimension(im,jm,kb) :: BOTRIB integer, dimension(im,jm) :: KPPbind REAL :: DEPTH, rhobot,ubot,vbot,WM,WS REAL :: DUZ, DVZ, RHOU, DRHOZ, RITOP, RIBOT REAL :: cff, cff_dn, cff_up, k_bl, dk_bl REAL :: GM1, dGM1DS,GT1,DGT1DS,A1,A2,A3,GM,GT,SIG,F1 integer :: indi,indj,indk,kbbl logical :: first, TRIG KPPbind(:,:)=kb kppb(:,:)=1.0 BOTRIB(:,:,:)=0.0 do indi=1,imm1 do indj=1,jmm1 ustarb2(indi,indj)= sqrt((.5e0*(wubot(indi,indj) $ * dum(indi,indj) + wubot(indi+1,indj) $ * dum(indi+1,indj)))**2 + $ (.5e0*(wvbot(indi,indj)*dvm(indi,indj) $ +wvbot(indi,indj+1)*dvm(indi,indj+1)))**2) hek(indi,indj) = 0.7 * sqrt(ustarb2(indi,indj)) $ / cor(indi,indj) enddo enddo do indj=1,JUP do indi=1,IUP first=.true. ! DEPTH = 0.0 - zz(kbm1) * h(indi,indj) ! call eqstate(S(indi,indj,kbm1),T(indi,indj,kbm1), $ 0.0,rhobot) ! if (dum(indi+1,indj).gt.0.1) then ubot = .5e0*(u(indi,indj,kbm1) $ +u(indi+1,indj,kbm1)) else ubot = .5e0*u(indi,indj,kbm1) endif ! if (dvm(indi,indj+1).gt.0.1) then vbot = .5e0*(v(indi,indj,kbm1) $ +v(indi,indj+1,kbm1)) else vbot = .5e0*v(indi,indj,kbm1) endif BOTRIB(indi,indj,:)=0.0 do indk = kbm2,2,-1 DEPTH = (zz(indk)-zz(kbm1)) * h(indi,indj) DUz = .5e0*(u(indi,indj,indk) $ +u(indi+1,indj,indk))-ubot DVz = .5e0*(v(indi,indj,indk) $ +v(indi,indj+1,indk))-vbot call eqstate(S(indi,indj,kbm1), $ T(indi,indj,kbm1),0.0,rhobot) call eqstate(S(indi,indj,indk), $ T(indi,indj,indk),0.0,RHOU) DRHOz=(RHOU)-rhobot RItop = -grav / rhoref * DRHOz * DEPTH RIbot= DUz**2 + DVz**2 BOTRIB(indi,indj,indk) = RItop / (RIbot + $ 1.e-8) !if (indi.eq.15 .and. indj.eq.15 .and. indk.eq.kbm2)then ! print*,'DU',DRHOz,DEPTH,RItop / (RIbot + 1.e-8) !endif enddo enddo enddo ! if (ionedim.eq.0) then do indk=1,kb call exchange2d_mpi(BOTRIB(:,:,indk),im,jm) enddo endif do indj=1,JUP do indi=1,IUP TRIG=.true. KPPb(indi,indj) = 1.0 if (botrib(indi,indj,kbm2).gt.KPP_RIc) then KPPb(indi,indj) = 0.0+ $ (KPP_Ric)*dzz(kbm2) $ *h(indi,indj)/(botRIB(indi,indj,kbm2)-0.0) if (h(indi,indj).lt.60..and.h(indi,indj).gt.20..and. $ botRIB(indi,indj,kbm2).lt.1.0) then ! print*,h(indi,indj),KPPb(indi,indj), ! $ botRIB(indi,indj,kbm2),botRIB(indi,indj,kb-3), ! $ dzz(kbm2),kpp_ric endif else do indk=kbm2,1,-1 if (BOTRIB(indi,indj,indk).gt.KPP_RIc .AND. trig) then KPPb(indi,indj) = (zz(indk)-zz(kbm1))* h(indi,indj) KPPb(indi,indj) = KPPb(indi,indj)- $ (BOTRIB(indi,indj,indk)-KPP_Ric)*dzz(indk) $ *h(indi,indj)/(botRIB(indi,indj,indk)- $ botRIB(indi,indj,indk-1)) trig=.false. elseif (zz(indk).gt.-0.5) then kppb(indi,indj)=h(indi,indj)*.5 trig=.false. endif enddo endif KPPb(indi,indj)=min(KPPb(indi,indj),h(indi,indj)) KPPb(indi,indj)=min(200.0,KPPb(indi,indj)) KPPbind(indi,indj)=kbm2 trig=.true. do indk=kb,2,-1 if (trig.and. $ abs((Z(indk)-Z(kb))*h(indi,indj)) $ .gt.KPPb(indi,indj)) then KPPbind(indi,indj)=indk trig=.false. endif enddo enddo enddo if (ionedim.eq.0) then Call exchange2d_mpi(KPPB,im,jm) call exchange2d_mpi(kppbind,im,jm) endif do indi=1,IUP do indj=1,JUP WM = kappa*sqrt(ustarb2(indi,indj)) WS = WM if (fsm(indi,indj).gt.0.1) then if (KPP_SIG_MATCH_BBL_LOG) then if (kppbind(indi,indj).lt.kbm1) then f1=0.0 kbbl=KPPbind(indi,indj) cff=1./abs(DZ(kbbl+1)*h(indi,indj)); cff_dn=(abs(Z(kbbl+1)*h(indi,indj)) $ -(h(indi,indj)-KPPb(indi,indj)))*cff cff_up=((h(indi,indj)-KPPb(indi,indj)) $ -abs(Z(kbbl)*h(indi,indj))) $ *cff K_bl = cff_dn*KMKPP(indi,indj,kbbl+1)+ $ cff_up*KMKPP(indi,indj,kbbl); dK_bl = cff*(KMKPP(indi,indj,kbbl+1) $ -KMKPP(indi,indj,kbbl)); Gm1 = K_bl/(abs(KPPb(indi,indj))*WM+small) dGm1dS = -dK_bl/(WM+small) + K_bl*f1; K_bl = cff_dn*KHKPP(indi,indj,kbbl+1) $ +cff_up*KHKPP(indi,indj,kbbl); dK_bl = cff*(KHKPP(indi,indj,kbbl+1) $ -KHKPP(indi,indj,kbbl)); Gt1 = K_bl/(abs(KPPb(indi,indj))*WS+small) dGt1dS = -dK_bl/(WS+small) + K_bl*f1; else GT1=0.0 GM1=0.0 endif endif do indk=1,kbm1 ! on z's (w levels) DEPTH = (z(indk) - z(kb)) * h(indi,indj) !is surf 0.0? et/el? IF (DEPTH.lt.KPPB(indi,indj)) then SIG=min(1.,DEPTH/KPPb(indi,indj)); A1 = SIG-2.0; A2 = 3.0-2.0*SIG; A3 = SIG-1.0; if (KPP_SIG_MATCH_BBL_LOG) THEN Gm = A1+A2*Gm1!+A3*dGm1dS; Gt = A1+A2*Gt1!+A3*dGt1dS; else GM=A1 GT=A1 endif kmkpp(indi,indj,indk)=DEPTH*WM*(1.+SIG*Gm) khkpp(indi,indj,indk)=DEPTH*WS*(1.+SIG*Gt) endif kmkpp(indi,indj,indk)=min(10.0, $ kmkpp(indi,indj,indk)) khkpp(indi,indj,indk)=min(10.0, $ khkpp(indi,indj,indk)) enddo else kmkpp(indi,indj,indk)=0.0 khkpp(indi,indj,indk)=0.0 endif enddo enddo if (ionedim.eq.0) then do indk=1,kb call exchange2d_mpi(KHKPP(:,:,indk),im,jm) call exchange2d_mpi(KMKPP(:,:,indk),im,jm) enddo endif end subroutine KPPbottom !/---------------------------------------------------------------------- !\---------------------------------------------------------------------- end module KPP