module physics use const use comm use var use io implicit real*8 (a-h,o-z) !---------------------------------------------------------------------- ! namelist !---------------------------------------------------------------------- logical, parameter :: lcfc = .false. ! flag for cfc module logical, parameter :: lkppmix = .true. ! flag for kpp vertical mixing scheme, parameters below are used logical, parameter :: lrich = .true. ! compute Ri-dependent mixing logical, parameter :: ltidal_mixing = .false. ! tidal mixing logical, parameter :: lhoriz_varying_bckgrnd = .false. ! horizontal varying backgrd vdc,vvc logical, parameter :: ldbl_diff = .false. ! dbl diffusion logical, parameter :: lcheckekmo = .false. ! check Ekman, Monin-Obhukov depth limit logical, parameter :: lshort_wave = .false. ! computing short-wave forcing logical, parameter :: linertial = .false. ! inertial mixing parameterization logical, parameter :: llangmuir = .false. ! Langmuir parameterization character*50 :: tidal_energy_file = 'data/tidal_enegry.nc' !---------------------------------------------------------------------- ! default coefficients and constants !---------------------------------------------------------------------- real*8, parameter :: convect_diff = 2000., &! diffusivity to mimic convection convect_visc = 2000. ! viscosity to mimic convection real*8, parameter :: epssfc = 0.1 ! non-dimensional extent of sfc layer !** parameters for Ri-dependent mixing integer, parameter :: num_v_smooth_Ri = 1 ! # of times to vertically smooth Ri real*8, parameter :: & Prandtl = 10.0, & rich_mix = 50.0, &! coefficient for rich number term Riinfty = 0.8, &! Rich. no. limit for shear instability BVSQcon = 0.0 ! Brunt-Vaisala square cutoff(s**-2) !** parameters for dbl diffusion real*8, parameter :: & Rrho0 = 2.55_r8, &! limit for double-diff density ratio dsfmax = 1.0_r8 ! max diffusivity for salt fingering !** parameters for bldepth real*8, parameter :: & vonkar = 0.4, &! von Karman constant cmonob = 1.0, &! coefficient for Monin-Obukhov depth cekman = 0.7, &! coefficient for Ekman depth concv = 1.7 ! minimum allowed value for buoyancy frequency ratio at entrainment depth real*8, parameter, dimension(k1) :: & Ricr = 0.3 ! crit Rich no. for bndy layer depth real*8, parameter, dimension(i0,j0) :: bolus_sp = 0.0 !scaled eddy-induced (bolus) speed used in inertial mixing !** parameters for velocity scale function (from Large et al.) real*8, parameter :: & zeta_m = -0.2, & zeta_s = -1.0, & c_m = 8.38, & c_s = 98.96, & a_m = 1.26, & a_s = -28.86 !** parameters for subroutine blmix: mixing within boundary layer real*8, parameter :: & cstar = 10.0 ! coeff for nonlocal transport !** parameters for tidal mixing real*8, parameter :: tidal_mix_max = 100. !** grid info for mixing calculation real*8,dimension(0:k0) :: zgrid,hwide !** background vvc and vdc real*8, dimension(i0,j0,k1) :: bckgrnd_vvc, bckgrnd_vdc !** mixing within boundary layer real*8 :: & cg, & !coefficient for counter-gradient term Vtc !resolution and buoyancy independent part of the !turbulent velocity shear coefficient (for bulk Ri no) !** Langmuir parameterization real*8, dimension(i0,j0) :: fstokes !ratio of stokes velocity to ustar !** coriolis coeff real*8, dimension(j0) :: fcort !** tidal mixing real*8, dimension(i0,j0,k1) :: tidal_diff real*8, dimension(i0,j0,k1) :: tidal_coef contains ! ---------------------------------------------------------------------- subroutine windglo(u,v,days,dt,odz) integer nld dimension u(i0,j0),v(i0,j0),odz(*) common/klimatglo/fnew,fold,new,nld common/windsglo/taux(i2,j2,12),tauy(i2,j2,12) ! the velocity arrays u,v have a perimeter "ghost zone". ! thus, we set wind forcing only in the interior zones. ! taux,tauy are surface wind stress components. ! taux,tauy units are force per unit area (i.e., energy per unit volume). ! tmp=odz(1)/rho ! all units are cgs, so we use rho=1. if(myid.eq.0) write(*,*) 'days = ',days tmp=dt*odz(1) do j=2,j1 do i=2,i1 u(i,j)=u(i,j)+tmp*(fold*taux(i-1,j-1,nld)+fnew*taux(i-1,j-1,new)) v(i,j)=v(i,j)+tmp*(fold*tauy(i-1,j-1,nld)+fnew*tauy(i-1,j-1,new)) enddo enddo end subroutine ! ---------------------------------------------------------------------- subroutine cfcsurf(c,t,s,days,dt,odz) dimension c(i0,j0),t(i0,j0),s(i0,j0),odz(*) !local logical,save :: first_call=.true. integer, parameter :: cfc_nt = 2 integer, parameter :: nyrbeg = 31 integer, parameter :: nyrend = 99 integer,parameter :: imt=i2 integer,parameter :: jmt=j2 real*8,save :: year(nyrend),p11(nyrend,cfc_nt),p12(nyrend,cfc_nt) real*8,save :: fice(i2,j2,12),xkw(i2,j2,12),patm(i2,j2,12) real*8,save :: ratio,day0,nday real*8,save :: p11_new(cfc_nt),p12_new(cfc_nt) real*8,save :: pnorth(cfc_nt),psouth(cfc_nt) real*8 :: rlon(imt,jmt),rlat(imt,jmt),CFCatm(imt,jmt,cfc_nt) real*8 :: Kw,CFCsat real*8 :: cfc_sc,cfc_alpha real*8, save :: CMAX,CMIN,STMAX,STMIN,SSMAX,SSMIN real*8, save :: KwMAX,scMAX,CFCMAX,alphaMAX,xkwMAX real*8, save :: KwMIN,scMIN,CFCMIN,alphaMIN,xkwMIN real*8 :: tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9 real*8 :: tmp10,tmp11,tmp12,tmp13,tmp14,tmp15,tmp16 INTEGER :: I,J,K,kn integer,save :: indx1,indx2,month !local if(first_call) then ! if(myid.eq.0)write(*,*)'days,dt=',days,dt if(myid.eq.0)write(*,*)'first call Read CFC gas and atm!' if(myid.eq.0)write(*,*)'i0,i1,i2=',i0,i1,i2 if(myid.eq.0)write(*,*)'j0,j1,j2=',j0,j1,j2 call read_gas(fice,xkw,patm) call read_cfcatm(cfc_nt,nyrbeg,nyrend,year,p11,p12) !Initial year in file indx1 = 90 !90 !31 !45 !1 !90 !we can use 29 =>spin up 2 yrs and run 37yrs indx2 = indx1+1 month=1 day0 = 0. !10800. if(myid.eq.0)write(*,*)'initial year,month:',1900+indx1,month STMAX = -9999.0 SSMAX = -9999.0 CMAX = -9999.0 KwMAX = -9999.0 scMAX = -9999.0 CFCMAX = -9999.0 alphaMAX = -9999.0 xkwMAX = -9999.0 STMIN = 9999.0 SSMIN = 9999.0 CMIN = 9999.0 KwMIN = 9999.0 scMIN = 9999.0 CFCMIN = 9999.0 alphaMIN = 9999.0 xkwMIN = 9999.0 do j=1,j2 do i=1,i2 rlon(i,j) = xdeg(i+1) rlat(i,j) = ydeg(j+1) ! if(myid.eq.0) write(99,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.1) write(98,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.2) write(97,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.3) write(96,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.4) write(95,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.5) write(94,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.6) write(93,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.7) write(92,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.8) write(91,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.9) write(90,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.10) write(89,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.11) write(88,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.12) write(87,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.13) write(86,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.14) write(85,*)rlon(i,j),rlat(i,j),xkw(i,j,12) ! if(myid.eq.15) write(84,*)rlon(i,j),rlat(i,j),xkw(i,j,12) enddo enddo first_call=.false. endif !first call nday = days-day0; if (mod(days,30.).eq.0.) then month = month +1 month = min(month,12) if(myid.eq.0)write(*,*)'month=',month end if if (mod(days,360.).eq.0.) then month = 1 day0 = days indx1 = indx2 indx2 = indx1+1 if(myid.eq.0)write(*,*)'New year: month,year1,year2=',month,indx1,indx2 end if ratio = nday/360. !360.!360.0 p11_new(1) = p11(indx1,1)*(1-ratio)+p11(indx2,1)*ratio p11_new(2) = p11(indx1,2)*(1-ratio)+p11(indx2,2)*ratio p12_new(1) = p11(indx1,1)*(1-ratio)+p11(indx2,1)*ratio p12_new(2) = p11(indx1,2)*(1-ratio)+p11(indx2,2)*ratio pnorth(1)=p11_new(1) pnorth(2)=p12_new(1) psouth(1)=p11_new(2) psouth(2)=p12_new(2) call cfc_interp(cfc_nt,pnorth, psouth, rlat, imt, jmt, CFCatm) ! if(mod(days,1.).eq.0) then ! do j=1,J2 ! do i=1,I2 ! if(myid.eq.0)write(99,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1) ! if(myid.eq.1)write(98,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1) ! if(myid.eq.2)write(97,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1) ! if(myid.eq.3)write(96,*)rlon(i,j),rlat(i,j),CFCatm(i,j,1) ! end do ! end do ! endif kn = 11 !11: CFC-11 12:CFC-12 ! if(myid.eq.0) write(*,*)'T,S=',T(70,60),S(70,60) ! if(myid.eq.0) call sc_cfc(T(70,60),kn,cfc_sc) ! if(myid.eq.0) call sol_cfc(T(70,60),S(70,60),kn,cfc_alpha) ! if(myid.eq.0)write(*,*)'alpha=',cfc_sc,cfc_alpha do j=1,j2 do i=1,i2 call sc_cfc(t(I+1,J+1),kn,cfc_sc) call sol_cfc(t(I+1,J+1),s(I+1,J+1),kn,cfc_alpha) CFCsat = cfc_alpha*CFCatm(I,J,1) !1=CFC-11 ! convert to pmol/kg for comparison with GLODAP data set CFCsat = CFCsat*1e12/1030. !pmol/kg Kw = (1-fice(I,J,month))*xkw(I,J,month)*(cfc_sc/660.0)**(-0.5) !m/s Kw = Kw*ODZ(1)*100. !the latter *100 because ODZ is in 1/cm rather than 1/m Kw = Kw*dt ! if(myid.eq.0) then ! if(i.eq.50.and.j.eq.50) then ! write(*,*)'check kw',myid,month,xkw(I+1,J+1,month) ! endif ! endif c(I+1,J+1) = c(I+1,J+1) + Kw*(CFCsat-c(I+1,J+1))*IN(I+1,J+1,1) STMAX=MAX(STMAX,T(I+1,J+1)) SSMAX=MAX(SSMAX,S(I+1,J+1)) CMAX=MAX(CMAX,C(I+1,J+1)) KwMAX=MAX(KwMAX,Kw) scMAX=MAX(scMAX,cfc_sc) CFCMAX=MAX(CFCMAX,CFCsat) alphaMAX=MAX(alphaMAX,cfc_alpha) xkwMAX=MAX(xkwMAX,xkw(I,J,month)) STMIN=MIN(STMIN,T(I+1,J+1)) SSMIN=MIN(SSMIN,S(I+1,J+1)) CMIN=MIN(CMIN,C(I+1,J+1)) KwMIN=MIN(KwMIN,Kw) scMIN=MIN(scMIN,cfc_sc) CFCMIN=MIN(CFCMIN,CFCsat) alphaMIN=MIN(alphaMIN,cfc_alpha) xkwMIN=MIN(xkwMIN,xkw(I,J,month)) end do end do tmp1=cmax tmp2=cmin tmp3=cfcmax tmp4=cfcmin tmp5=kwmax tmp6=kwmin tmp7=scmax tmp8=scmin tmp9=alphamax tmp10=alphamin tmp11=stmax tmp12=stmin tmp13=ssmax tmp14=ssmin tmp15=xkwmax tmp16=xkwmin CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp1,cmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp2,cmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp3,cfcmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp4,cfcmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp5,kwmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp6,kwmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp7,scmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp8,scmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp9,alphamax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp10,alphamin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp11,stmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp12,stmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp13,ssmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp14,ssmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp15,xkwmax,1,MPI_REAL8,MPI_MAX,M_CART,IERR) CALL MPI_BARRIER(M_CART,IERR) CALL MPI_ALLREDUCE(tmp16,xkwmin,1,MPI_REAL8,MPI_MIN,M_CART,IERR) if(myid.eq.0) then if (mod(days,1.0).eq.0.) then write(*,*)'CFC check! DT=',DT write(*,*)days,month,nday,day0,ratio,indx1,indx2 write(*,*)stmax,ssmax,cmax,cfcmax,kwmax,scmax,alphamax write(*,*)stmin,ssmin,cmin,cfcmin,kwmin,scmin,alphamin write(*,*)xkwmax,xkwmin,odz(1),odz(1)*100.0*dt end if endif end subroutine subroutine sc_cfc(t,kn,cfc_sc) !--------------------------------------------------- ! CFC 11 and 12 schmidt number ! as a fonction of temperature. ! ! ref: Zheng et al (1998), JGR, vol 103,No C1 ! ! t: temperature (degree Celcius) ! kn: = 11 for CFC-11, 12 for CFC-12 ! ! J-C Dutay - LSCE !--------------------------------------------------- ! use cfc_glbl, only: rkind implicit none integer, intent(in) :: kn real*8 :: a1( 11: 12), a2( 11: 12), a3( 11: 12), a4( 11: 12) real*8,intent(in) :: t real*8,intent(out) :: cfc_sc ! ! coefficients with t in degre Celcius ! ------------------------------------ a1(11) = 3501.8 a2(11) = -210.31 a3(11) = 6.1851 a4(11) = -0.07513 ! a1(12) = 3845.4 a2(12) = -228.95 a3(12) = 6.1908 a4(12) = -0.067430 ! cfc_sc = a1(kn)+a2(kn)*t+a3(kn)*t*t+a4(kn)*t*t*t end subroutine subroutine sol_cfc(pt,ps,kn,cfc_alpha) implicit none real*8,intent(in) :: pt, ps real*8,intent(out) :: cfc_alpha real*8 :: ta,d real*8 :: a1 ( 11: 12), a2 ( 11: 12), a3 ( 11: 12), a4 ( 11: 12) real*8 :: b1 ( 11: 12), b2 ( 11: 12), b3 ( 11: 12) integer, intent(in) :: kn !! !! coefficient for solubility in mol/l/atm !! ---------------------------------------- ! ! for CFC 11 ! ---------- a1 ( 11) = -229.9261 a2 ( 11) = 319.6552 a3 ( 11) = 119.4471 a4 ( 11) = -1.39165 b1 ( 11) = -0.142382 b2 ( 11) = 0.091459 b3 ( 11) = -0.0157274 ! ! for CFC/12 ! ---------- a1 ( 12) = -218.0971 a2 ( 12) = 298.9702 a3 ( 12) = 113.8049 a4 ( 12) = -1.39165 b1 ( 12) = -0.143566 b2 ( 12) = 0.091015 b3 ( 12) = -0.0153924 ! ta = ( pt + 273.16)* 0.01 d = ( b3 ( kn)* ta + b2 ( kn))* ta + b1 ( kn) ! ! cfc_alpha = exp ( a1 ( kn) + a2 ( kn)/ ta + a3 ( kn)* log ( ta ) + & a4 (kn)* ta * ta + ps* d ) ! ! conversion from mol/(l * atm) to mol/(m^3 * atm) ! ------------------------------------------------ cfc_alpha = 1000. * cfc_alpha ! ! conversion from mol/(m^3 * atm) to mol/(m3 * pptv) ! -------------------------------------------------- cfc_alpha = 1.0e-12 * cfc_alpha end subroutine ! ---------------------------------------------------------------------- subroutine cfc_interp(cfc_nt,pnorth,psouth,rlat,imt,jmt,CFCatm) integer, parameter :: ijmax=200*200 integer, intent(in) :: cfc_nt, imt, jmt integer,save :: ientry integer :: i, j, ij, n real*8,save :: ys, yn real*8,save :: xintp(ijmax) real*8, intent(in) :: Pnorth(cfc_nt), Psouth(cfc_nt) real*8, intent(in) :: rlat(imt,jmt) real*8, intent(out) :: CFCatm(imt,jmt,cfc_nt) DATA ys, yn / -10., 10./ DATA ientry /0/ ientry = ientry + 1 ! Test to see if ijmax is large enough if (ijmax .lt. imt*jmt) then print *," cfc_interp: ERROR -> ijmax must be at least",imt*jmt stop endif ! IF Block to be executed ONLY on first call to this routine if (ientry .eq. 1) then do j = 1,jmt do i = 1,imt ij = (j-1)*imt + i if (rlat(i,j) .ge. yn) then xintp(ij) = 1. else if (rlat(i,j) .le. ys) then xintp(ij) = 0. else xintp(ij) = (rlat(i,j) - ys)/(yn-ys) endif end do end do endif ! Block to be executed every pass do n=1,cfc_nt do j=1,jmt do i=1,imt ij = (j-1)*imt + i CFCatm(i,j,n) = xintp(ij) * Pnorth(n) + & (1.0 - xintp(ij)) * Psouth(n) end do end do end do end subroutine ! ---------------------------------------------------------------------- subroutine read_cfcatm(cfc_nt,nyrbeg,nyrend,year,p11,p12) ! needs modification to dbl !local integer :: i, n integer, parameter :: iu = 10 integer,intent(in) :: nyrbeg,nyrend,cfc_nt real*8,intent(out) :: year(nyrend), p11(nyrend,cfc_nt), p12(nyrend,cfc_nt) character(len=80) :: filen year = 0. p11 = 0. p12 = 0. ! ! OPEN FILE ! --------- filen='data/cfc1112_new.atm' open(UNIT=iu,FILE=filen,STATUS='old') ! ! Skip over 1st six descriptor lines ! ---------------------------------- do i=1,6 read(iu,99) end do ! ! READ FILE ! --------- do n = nyrbeg, nyrend READ(iu,*)year(n), p11(n,1), p12(n,1), p11(n,2), p12(n,2) year(n) = year(n) - 0.5 if(myid.eq.0) then write(*,*) 'read cfcatm file', year(n),p11(n,1), p12(n,1), p11(n,2), p12(n,2) endif end do close(iu) 99 FORMAT(1x) 100 FORMAT(f7.2, 4f8.2) end subroutine ! ---------------------------------------------------------------------- subroutine pp82glo ! ---------------------------------------------------------------------- ! vertical mixing coefficients ! approach: ! set background vertical diffusivities to dmz0 (o(0.1) cm-cm/sec) and ! add richardson number based vertical diffusivity plus numerical ! contribution according to vertical cell reynolds number. ! vertical cell re is o(100), except during winter cooling conditions ! and in wind-blown surface mixed layer, because internal waves, which ! dominate w below the sml in the model results during summer, do not ! mix t or s in nature. ! note: one cannot depend on diffusive closure by itself if one wants to ! model possible contra-diffusive effe! parameter(i01=i0+1) common/windmxglo/vbk(k2),hbk(k2) common/tauglo/add(i2,j2,k2) do k=1,k2 l=k+1 dzwl=orzmx/odzw(l) ! stability limit emax=.4/(dt*odzw(l)**2) do j=2,j1 do i=2,i1 ! pacanowski and philander (1981): evisc=a*(1/(1+b*ri))**n ! evisc=turbulent eddy viscosity; diffuse=turbulent diffusivity ! ri=max(0.,980.*(rho(i,j,l)-rho(i,j,k))*odzw(l)/(odzw(l)**2* ri=max(-0.02d0,980.*(rho(i,j,l)-rho(i,j,k))*odzw(l)/(odzw(l)**2*(0.001+(u2(i,j,l)-u2(i,j,k))**2+(v2(i,j,l)-v2(i,j,k))**2))) rrinv=1./(1.+5.*ri) tmp=rrinv**2 ! ri- and cell re- dependent vertical mixing ! numerically kosher because, for large ri, laminar flow dominates ! and time mean overshoot effe! are limited by time mean w being small temp=tmp*abs(w(i,j,l))*dzwl evisc=5.*tmp diffuse=evisc*rrinv ! add parameterizes extra momentum dissipation along steep slopes, ! high and low latitudes, and equatorial regions due to breaking of ! waves that were generated by big storm events that are missing when ! using climatological winds. add includes vbk(k) term. ev(i-1,j-1,k)=evisc+add(i-1,j-1,k)+temp hv(i-1,j-1,k)=diffuse+hbk(k)+temp ! apply stability limit and odzw factor tmp=odzw(l)*iw(i,j,l) ev(i-1,j-1,k)=tmp*min(emax,ev(i-1,j-1,k)) hv(i-1,j-1,k)=tmp*min(emax,hv(i-1,j-1,k)) enddo enddo enddo end subroutine ! ---------------------------------------------------------------------- subroutine pp82glonew ! ---------------------------------------------------------------------- ! vertical mixing coefficients ! approach: ! set background vertical diffusivities to dmz0 (o(0.1) cm-cm/sec) and ! add richardson number based vertical diffusivity plus numerical ! contribution according to vertical cell reynolds number. ! vertical cell re is o(100), except during winter cooling conditions ! and in wind-blown surface mixed layer, because internal waves, which ! dominate w below the sml in the model results during summer, do not ! mix t or s in nature. ! note: one cannot depend on diffusive closure by itself if one wants to ! model possible contra-diffusive effe! parameter(i01=i0+1) common/windmxglo/vbk(k2),hbk(k2) rate=-.05 jt=mylat*j2 do k=1,k2 l=k+1 dzwl=orzmx/odzw(l) ! deepen and augment mixing in high lats to emulate big storms tmpk=10.*exp(-.0001*z(2*k+1)) ! stability limit emax=.4/(dt*odzw(l)**2) elim=min(200.,emax) do j=2,j1 ! decrease background mixing except in high lats to emulate big ! storms within 20 grid intervals (~2 deg lat) of lat boundaries tmp=tmpk*(exp(rate*abs(2-j-jt)))+exp(rate*abs(j1t-j-jt)) vb=vbk(k)+tmp hb=hbk(k)+tmp do i=2,i1 ! pacanowski and philander (1981): evisc=a*(1/(1+b*ri))**n ! evisc=turbulent eddy viscosity; diffuse=turbulent diffusivity tmp=odzw(l) ! starting at day 0 of yr 16, decrease ri-based vertical mixing ! gap in izu ridge is also deepened to emulate unresolved deep kuroshio paths ! thus allowing it to pass over the ridgeline with less anomalous vortex ! squashing and thus taking a non-meander path when appropriate. oldfact=980.*(rho(i,j,l)-rho(i,j,k))*tmp/(tmp**2*(0.001+(u2(i,j,l)-u2(i,j,k))**2+(v2(i,j,l)-v2(i,j,k))**2)) ! big ri decreases vertical mixing ri=max(-0.19,5.*oldfact) rrinv=1./(1.+5.*ri) tmp=rrinv**2 ! ri- and cell re- dependent vertical mixing ! numerically kosher because, for large ri, laminar flow dominates ! and time mean overshoot effe! are limited by time mean w being small temp=.01*tmp*abs(w(i,j,l))*dzwl evisc=tmp diffuse=evisc*rrinv ! evisc, diffuse, hb and vb units are cm-cm/sec ev(i-1,j-1,k)=evisc+vb+temp hv(i-1,j-1,k)=diffuse+hb+0.1*temp ! apply stability limit and odzw factor tmp=odzw(l)*iw(i,j,l) ev(i-1,j-1,k)=tmp*min(elim,ev(i-1,j-1,k)) hv(i-1,j-1,k)=tmp*min(elim,hv(i-1,j-1,k)) enddo enddo enddo ! ====================================================================== ! add vertical mixing to emulate effect of equatorial vortices and ! penetrative convection in acc. ! ====================================================================== ! stacked vortices in 10-point equator band (hardwired to 1/4 deg res) ! strong velocity oscillations having ~20-day time scale are observed ! within a few deg of equator ! 33-point equatorial and acc bands ! ! do 770 i=1,i2 ! do 770 j=1,j2 ! do 770 k=1,k2 ! 770 scr(i+1,j+1,k)=ev(i,j,k) sc1=.01 sc2=.01 do k=1,k2 l=k+1 ! 100m scale height of added vertical mixing tmp=exp(-.0001*z(2*k+1)) do j=1,j2 tmp1=10./(1.+sc1*((0-ydeg(j))/(ydeg(j+1)-ydeg(j)))**2) addd=tmp1*odzw(l)*tmp do i=1,i2 ! insulated bottom, nonlinear drag invoked in fs ad=addd*iw(i+1,j+1,l) hv(i,j,k)=hv(i,j,k)+ad ev(i,j,k)=ev(i,j,k)+ad enddo enddo enddo end subroutine ! ---------------------------------------------------------------------- subroutine smagglo ! ---------------------------------------------------------------------- ! horizontal mixing coefficients ! approach:smagorinsky, j. general circulation experiments with the ! primitive equations, i. the basi! experiment. monthly weather rev. 1963, ! 91, 99-164. parameter(i01=i0+1) do k=1,k1 do j=2,j1 temp=ocs(j)*ody(j) do i=2,i1 scr(i,j,k)=0.12*dx(j)*dy(j)*in(i,j,k)*sqrt(((u(i,j,k)-u(i-1,j,k))*odx(j))**2+((csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)**2) enddo enddo enddo call mpi_sendrecv(scr(2,1,1),1,m_vlon,m_w,1,scr(i0,1,1),1,m_vlon,m_e,1,m_cart,istat,ierr) call mpi_sendrecv(scr(i1,1,1),1,m_vlon,m_e,1,scr(1,1,1),1,m_vlon,m_w,1,m_cart,istat,ierr) call mpi_sendrecv(scr(1,2,1),1,m_vlat,m_s,1,scr(1,j0,1),1,m_vlat,m_n,1,m_cart,istat,ierr) call mpi_sendrecv(scr(1,j1,1),1,m_vlat,m_n,1,scr(1,1,1),1,m_vlat,m_s,1,m_cart,istat,ierr) do k=1,k1 do j=1,j0 do i=1,i1 dmx(i,j,k)=dm0+.5*(scr(i,j,k)+scr(i+1,j,k)) enddo enddo do i=1,i0 dmy(i,1,k)=dm0+scr(i,2,k) dmy(i,j1,k)=dm0+scr(i,j2,k) do j=1,jf ! note: except in hi lats, dmy slightly violates ang. mom. conservation dmy(i,j,k)=dm0+.5*(scr(i,j,k)+scr(i,j+1,k)) enddo enddo enddo end subroutine ! ---------------------------------------------------------------------- subroutine salsurfglo(s) real*8, dimension(i0,j0), intent(inout) :: s real*8 dtmp ! Salinity= Salinity + [ tmp * sdot ] (sdot==evapo) ! tmp converts precipitation and evaporation rate (kg/m2/s) to psu ! hence tmp is DT/DZ /waterdensity * salinity ! tmp = dt*odz/1000 * psu dtmp=dt*odz(1)/1000.d0*100.d0 !the latter *100 because ODZ is in 1/cm rather than 1/m DO j=1,j2 DO i=1,i2 ! S(I+1,J+1)=S(I+1,J+1)*(1-(DTMP*(EVAPO(i,j,1)+RAIN(I,J,1)))) + ! 1 DTMP*SALT(i,j,1) s(i+1,j+1)=s(i+1,j+1)-34.7d0*(dtmp*(evapo(i,j)+prec(i,j))) & + dtmp*salt(i,j) ! EVAPO<0 (negative up); PREC>0 (positive down); - (negative EVAPO+PREC) --> ! SALINITY INCREASE ! unit of SALT is kg(psu)/m^2/s enddo enddo end subroutine salsurfglo ! ---------------------------------------------------------------------- ! ---------------------------------------------------------------------- subroutine qsurfglo(t) real*8, dimension(i0,j0), intent(inout) :: t real*8 :: dtmp ! the temperature array t has a permiter "ghost zone". ! thus, we set atmospheri! heat exchange only in the interior zones. ! typical latitudinal range: -50 to 50 watts per square meter. ! 1 watt per square meter = 1 erg per second per equare cm. ! example: surface heating qdot = 50 watts per square meter !qdot=50. ! tmp is conversion factor from watts per square meter ! to deg ! per time step in top model layer. ! tmp=dt*odz(1)/rho/cp ! all units are cgs, so we use rho=1 !dog -> 1.026 ! cp for water is 2.5e4 ergs/gram/deg c. ! Cp for water is 3.996 !POP_CpSW = 3.996e7_POP_r8 ! specific heat salt water !POP_CpSW = SHR_CONST_CPSW*10000.0_POP_r8 ! erg/g/K dtmp=dt*odz(1)/(3.996d4*1.026) do j=1,j2 do i=1,i2 ! add because heat is going IN (positive downward) ! t(i+1,j+1)=t(i+1,j+1)+dtmp*(qdot(i,j,1)+qdot2(i,j,1)) t(i+1,j+1)=t(i+1,j+1)+dtmp*(qdot(i,j)) enddo enddo end subroutine qsurfglo ! ---------------------------------------------------------------------- !subroutine sunglo(t,dt,odz,qprof,kb) ! integer*2 kb ! dimension t(i0,j0,k1),odz(k1),qprof(k1),kb(i0,j0) ! this is like subroutine qsurf except ! we use predetermined input array qprof to vertically distribute qdot. ! qprof(k) is the cumulative fraction of total heat absorbed to depth ! of bottom of layer k for a pure water column. ! example: surface radiation flux of 10 watts per square meter ! qdot=10. ! fract=qprof(1) ! do k=1,k1 ! dtemp=fract*qdot*dt/2.5e4 ! odz(k) is not the correct factor for layers deeper than the local bottom. ! thus, odz factor must be added inside loop 99. ! this quickly programmed logi! must be checked. ! do j=1,j2 ! do i=1,i2 ! kk=min(kb(i+1,j+1),k) ! t(i+1,j+1,kk)=t(i+1,j+1,kk)+dtemp*odz(kk) ! enddo ! enddo ! if (k.lt.k1) fract=qprof(k+1)-fract ! enddo !end subroutine !-------------------------------------------------------------------------- subroutine kppsrcglo(t,s) real*8 :: dtmp1, dtmp2 real*8, dimension(i0,j0,k1), intent(inout) :: t,s ! KPP_SRC unit: T (degree C/s) S (g/kg/S) ! SO ONLY NEED DT CONVERSION ! NEED TO CHECK THE UNIT OF STF/DZ*VDC again dtmp1 = dt dtmp2 = dt do k=1,k1 do j=1,j2 do i=1,i2 t(i+1,j+1,k)=t(i+1,j+1,k) + dtmp1*kpp_src(i+1,j+1,k,1)*IN(i+1,j+1,k) s(i+1,j+1,k)=s(i+1,j+1,k) + dtmp2*kpp_src(i+1,j+1,k,2)*IN(i+1,j+1,k) !add the KPP NON-LOCAL SOURCE enddo enddo enddo end subroutine kppsrcglo subroutine kppglo implicit none integer :: i,j,k,l real*8 :: emax real*8, dimension(i0,j0,k1) :: vvc real*8, dimension(i0,j0,0:k0,2) :: vdc ! !DESCRIPTION: ! This is the main driver routine which calculates EV,HV and KPP_SRC ! for the KPP mixing scheme as outlined in ! Large, McWilliams and Doney, Reviews of Geophysics, 32, 363 ! (November 1994). KPP_SRC part is treated in ==surface effect== do j = 1,j2 do i = 1,i2 STF(i+1,j+1,1)=IN(i+1,j+1,1)*qdot(i,j)/(3.996d4*1.026) ! W/m^2 to deg C*cm/s STF(i+1,j+1,2)=IN(i+1,j+1,1)*((evapo(i,j)+prec(i,j))*-34.7d0*1.e-4_r8+ & ! salt(i,j)*1.e-3_r8*0.1_r8) !kg/m^2/s to psu*cm/s SHF_QSW(i+1,j+1)=IN(i+1,j+1,1)*qdot2(i,j)/(3.996d4*1.026) ! W/m^2 to deg C*cm/s SMFT(i+1,j+1,1)=IN(i+1,j+1,1)* ( fold*taux(i,j,nld) + fnew*taux(i,j,new) ) !dyne/cm^2 SMFT(i+1,j+1,2)=IN(i+1,j+1,1)* ( fold*tauy(i,j,nld) + fnew*tauy(i,j,new) ) !dync/cm^2 enddo enddo STF(:,:,2) = 0. CALL MPI_EXCH_2D_R8(STF(:,:,1)) CALL MPI_EXCH_2D_R8(STF(:,:,2)) CALL MPI_EXCH_2D_R8(SHF_QSW(:,:)) CALL MPI_EXCH_2D_R8(SMFT(:,:,1)) CALL MPI_EXCH_2D_R8(SMFT(:,:,2)) call qvolumeglo(t2,trans,0) !! flag = 0 ,calculate trans only call mpi_exch_3d_r8(u2) call mpi_exch_3d_r8(v2) call mpi_exch_3d_r8(t2) call mpi_exch_3d_r8(s2) trcr(:,:,:,1)=t2 trcr(:,:,:,2)=s2*1.e-3_r8 ! change to pop2 salt unit (g/g) call vmix_coeffs_kpp(IN,trans,kpp_src, vdc, vvc, trcr, u2,v2,rho,& stf, shf_qsw, convect_diff, convect_visc, & smft=smft) ! add = 0. do k=1,k2 l=k+1 emax=.5d0/(dt*odzw(l)**2) emax=MIN(1200.d0,emax) do j=1,j2 do i=1,i2 dtmp=odzw(l)*iw(i+1,j+1,l) !! limit the value not to over 900 or .5/(DT*ODZW(L)**2) ev(i,j,k)=dtmp*MAX(MIN(emax,vvc(i+1,j+1,k)+add(i,j,k)),vbk(k)*5.d0) hv(i,j,k)=dtmp*MAX(MIN(emax,vdc(i+1,j+1,k,1)+add(i,j,k)),hbk(k)*5.d0) end do end do end do end subroutine kppglo !-------------------------------------------------------------------------- subroutine vmix_coeffs_kpp(IN,trans,kpp_src, vdc, vvc, trcr, uuu,vvv,rhomix,& stf, shf_qsw, convect_diff, convect_visc, & smf,smft) ! !DESCRIPTION: ! This is the main driver routine which calculates the vertical ! mixing coefficients for the KPP mixing scheme as outlined in ! Large, McWilliams and Doney, Reviews of Geophysics, 32, 363 ! (November 1994). The non-local mixing is also computed here, but ! is treated as a source term in baroclinic. implicit none ! !INPUT PARAMETERS: integer*2, dimension(i0,j0,k1), intent(in) :: IN real*8, dimension(i0,j0,k1,2), intent(in) :: & trcr ! tracers at current time real*8, dimension(i0,j0,k1), intent(in) :: & uuu,vvv, & ! velocities at current time trans ! transmission (1 to 0) from surface (calculated in kppglo) real*8, dimension(i0,j0,k1), intent(in) :: & rhomix ! density at mix time real*8, dimension(i0,j0,2), intent(in) :: & stf ! surface forcing for all tracers real*8, dimension(i0,j0), intent(in) :: & shf_qsw ! short-wave forcing real*8, dimension(i0,j0,2), intent(in), optional :: & smf, &! surface momentum forcing at U points smft ! surface momentum forcing at T points ! *** either one or the other (not ! *** both) should be passed real*8, intent(in) :: & convect_diff, &! diffusivity to mimic convection convect_visc ! viscosity to mimic convection ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0,k1), intent(inout) :: & vvc ! viscosity for momentum diffusion real*8, dimension(i0,j0,0:k0,2),intent(inout) :: & vdc ! diffusivity for tracer diffusion real*8, dimension(i0,j0,k1,2), intent(out) :: & kpp_src !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & k, &! vertical level index i,j, &! horizontal loop indices n, &! tracer index mt2 ! index for separating temp from other trcrs integer, dimension(i0,j0) :: & kbl ! index of first lvl below hbl real*8, dimension(i0,j0) :: & ustar, &! surface friction velocity bfsfc, &! surface buoyancy forcing work1,work2,&! temporary storage fcon, &! convection temporary stable ! = 1 for stable forcing; = 0 for unstable forcing real*8, dimension(i0,j0,k1) :: & dbloc, &! buoyancy difference between adjacent levels dbsfc, &! buoyancy difference between level and surface ghat ! non-local mixing coefficient real*8, dimension(i0,j0,0:k0) :: & visc ! local temp for viscosity !dog !----------------------------------------------------------------------- ! initialize !----------------------------------------------------------------------- call init_vmix_kpp !----------------------------------------------------------------------- ! compute buoyancy differences at each vertical level. ! derived variables: dbloc, dbsfc !----------------------------------------------------------------------- call buoydiff(dbloc, dbsfc, trcr) !----------------------------------------------------------------------- ! compute mixing due to shear instability, internal waves and ! convection ! derived variables: vdc,visc !----------------------------------------------------------------------- call ri_iwmix(dbloc, visc, vdc, uuu, vvv) !----------------------------------------------------------------------- ! compute double diffusion if desired !----------------------------------------------------------------------- if (ldbl_diff) call ddmix(vdc,trcr) !----------------------------------------------------------------------- ! compute boundary layer depth ! derived variables: kpp_hblt, ustar, bfsfc, stable, kbl !----------------------------------------------------------------------- if (present(smft)) then call bldepth (trans,dbloc, dbsfc, trcr, uuu, vvv, stf, shf_qsw, & kpp_hblt, ustar, bfsfc, stable, kbl, & smft=smft) else call bldepth (trans,dbloc, dbsfc, trcr, uuu, vvv, stf, shf_qsw, & kpp_hblt, ustar, bfsfc, stable, kbl, & smf=smf) endif !----------------------------------------------------------------------- ! compute boundary layer diffusivities ! derived variables: visc,vdc,ghat !----------------------------------------------------------------------- call blmix(visc, vdc, kpp_hblt, ustar, bfsfc, stable, & kbl, ghat) !----------------------------------------------------------------------- ! ! consider interior convection: ! ! compute function of Brunt-Vaisala squared for convection. ! ! use either a smooth ! ! WORK1 = N**2, FCON is function of N**2 ! FCON = 0 for N**2 > 0 ! FCON = [1-(1-WORK1/BVSQcon)**2]**3 for BVSQcon < N**2 < 0 ! FCON = 1 for N**2 < BVSQcon ! ! or a step function. The smooth function has been used with ! BVSQcon = -0.2e-4_dbl_kind. ! ! after convection, average viscous coeffs to U-grid and reset sea ! floor values ! !----------------------------------------------------------------------- do k=1,k1-1 work1 = dbloc(:,:,k)/(zgrid(k) - zgrid(k+1)) if (BVSQcon /= 0.0) then work2 = min(1.0_r8-(max(work1,BVSQcon))/BVSQcon, 1.0_r8) fcon = (1.0_r8 - work2*work2)**3 else where (work1 > 0.0_r8) fcon = 0.0_r8 elsewhere fcon = 1.0_r8 end where endif !*** add convection and reset sea floor values to zero do j=1,j0 do i=1,i0 if ( k >= kbl(i,j) ) then visc(i,j,k) = visc(i,j,k) + convect_visc * fcon(i,j) vdc(i,j,k,1) = vdc(i,j,k,1) + convect_diff * fcon(i,j) vdc(i,j,k,2) = vdc(i,j,k,2) + convect_diff * fcon(i,j) endif if ( k >= int(kb(i,j)) ) then visc(i,j,k ) = 0.0_r8 vdc (i,j,k,1) = 0.0_r8 vdc (i,j,k,2) = 0.0_r8 endif end do end do vvc(:,:,k) = merge(visc(:,:,k), 0.0_r8, ( k < int(kb(:,:)) ) ) enddo VDC(:,:,k1,:) = 0. VVC(:,:,k1) = 0. !----------------------------------------------------------------------- ! add ghatp term from previous computation to right-hand-side ! source term on current row !----------------------------------------------------------------------- do n=1,2 mt2=min(n,2) kpp_src(:,:,1,n) = stf(:,:,n)*odz(1) & *(-vdc(:,:,1,mt2)*ghat(:,:,1)) do k=2,k1 kpp_src(:,:,k,n) = stf(:,:,n)*odz(k) & *( vdc(:,:,k-1,mt2)*ghat(:,:,k-1) & -vdc(:,:,k ,mt2)*ghat(:,:,k )) enddo enddo end subroutine vmix_coeffs_kpp !-------------------------------------------------------------------------- subroutine qvolumeglo(t,transmitted,flag) ! tmp is conversion factor from watts per square meter ! to deg ! per time step in top model layer. ! tmp=dt*odz(1)/rho/cp ! all units are cgs, so we use rho=1 !dog -> 1.026 ! cp for water is 2.5e4 ergs/gram/deg c. ! Cp for water is 3.996 implicit none real*8, dimension(i0,j0,k1), intent(inout) :: t,transmitted integer, intent(in) :: flag !-------------------------------------------------------------------------- ! local variables !-------------------------------------------------------------------------- real*8 :: myZ,R,xi1,xi2,total_trans,dtmp real*8 :: qdot2sum,qdot2ori,qdotori integer :: i,j,k R=.58d0 xi1=.35d0 xi2=23.d0 do j=1,j0 do i=1,i0 total_trans=0.d0 do k=1,k1 myZ=-z(2*k)/100.d0 transmitted(i,j,k)=(R*exp(myZ/xi1)+(1-R)*exp(myZ/xi2))* & in(i,j,k)*odz(k) total_trans=total_trans+transmitted(i,j,k) enddo do k=1,k1 if (total_trans.gt.1.d-6) then transmitted(i,j,k)=transmitted(i,j,k)/total_trans else transmitted(i,j,k)=0.d0 endif enddo enddo enddo if (flag==0) goto 8888 ! absorption*QDOT is added to T qdot2ori=SUM(qdot2(:,:)*in(2:i1,2:j1,1)) dtmp=qdot2ori CALL MPI_ALLREDUCE(DTMP,QDOT2ORI,1,MPI_REAL8,MPI_SUM,M_CART,IERR) qdotori=SUM(qdot(:,:)*in(2:i1,2:j1,1)) dtmp=qdotori CALL MPI_ALLREDUCE(DTMP,QDOTORI,1,MPI_REAL8,MPI_SUM,M_CART,IERR) qdot2sum=0.d0 do k=1,K1 dtmp=dt*odz(k)/(3.996d4*1.026) do j=1,j2 do i=1,i2 t2(i+1,j+1,k)=t2(i+1,j+1,k) + & dtmp*qdot2(i,j)*transmitted(i+1,j+1,k)*IN(i+1,j+1,k) !add because heat is going IN qdot2sum=qdot2sum+qdot2(i,j)*transmitted(i+1,j+1,k)*IN(i+1,j+1,k) enddo enddo enddo dtmp=qdot2sum CALL MPI_ALLREDUCE(DTMP,QDOT2SUM,1,MPI_REAL8,MPI_SUM,M_CART,IERR) if(myid==0)WRITE(*,*)'TOTALQ:',QDOTORI,QDOT2ORI,QDOT2SUM 8888 continue end subroutine qvolumeglo !-------------------------------------------------------------------------- !-------------------------------------------------------------- subroutine buoydiff(dbloc,dbsfc,trcr) ! !DESCRIPTION: ! This routine calculates the buoyancy differences at model levels. ! ! !REVISION HISTORY: ! same as module ! implicit none ! !INPUT PARAMETERS: real*8, dimension(i0,j0,k1,2), intent(in) :: & trcr ! tracers at current time ! !OUTPUT PARAMETERS: real*8 , dimension(i0,j0,k1), intent(out) :: & dbloc, &! buoyancy difference between adjacent levels dbsfc ! buoyancy difference between level and surface !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & k, &! vertical level index i,j, &! horizontal indices kprev, klvl, ktmp ! indices for 2-level TEMPK array integer, dimension(i0,j0) :: kmt real*8 , dimension(i0,j0) :: & rho1, &! density of sfc t,s displaced to k rhokm, &! density of t(k-1),s(k-1) displaced to k rhok, &! density at level k tempsfc ! adjusted temperature at surface real*8 , dimension(i0,j0,2) :: & tempk ! temp adjusted for freeze at levels k,k-1 !----------------------------------------------------------------------- ! calculate density and buoyancy differences at surface !----------------------------------------------------------------------- tempsfc = merge(-2._r8,trcr(:,:,1,1),trcr(:,:,1,1) < -2._r8) klvl = 2 kprev = 1 tempk(:,:,kprev) = tempsfc dbsfc(:,:,1) = 0.0_r8 !----------------------------------------------------------------------- ! calculate DBLOC and DBSFC for all other levels !-----------------------------------------------------------------------\ do k = 2,k1 tempk(:,:,klvl) = merge(-2.0_r8,trcr(:,:,k,1),trcr(:,:,k,1) < -2.0_r8) ! calculate the density from Temperature and Salinity using an equation of state ! derived from McDougall, Wright, Jackett and Feistel (hereafter MWJF, 2001 ) call state(k, k, tempsfc, trcr(:,:,1,2), & rhofull=rho1) call state(k, k, tempk(:,:,kprev), trcr(:,:,k-1,2), & rhofull=rhokm) call state(k, k, tempk(:,:,klvl), trcr(:,:,k ,2), & rhofull=rhok) do j=1,j0 do i=1,i0 if (rhok(i,j) /= 0.0_r8) then dbsfc(i,j,k) = grav*(1.0_r8 - rho1 (i,j)/rhok(i,j)) dbloc(i,j,k-1) = grav*(1.0_r8 - rhokm(i,j)/rhok(i,j)) else dbsfc(i,j,k) = 0.0_r8 dbloc(i,j,k-1) = 0.0_r8 endif if ( k-1 >= int(kb(i,j)) ) dbloc(i,j,k-1) = 0.0_r8 enddo enddo ktmp = klvl klvl = kprev kprev = ktmp enddo dbloc(:,:,k1) = 0.0_r8 end subroutine buoydiff !-------------------------------------------------------------- !-------------------------------------------------------------- subroutine state(k, kk, tempk, saltk, & rhoout, rhofull, drhodt, drhods) ! !DESCRIPTION: ! Returns the density of water at level k from equation of state ! $\rho = \rho(d,\theta,S)$ where $d$ is depth, $\theta$ is ! potential temperature, and $S$ is salinity. the density can be ! returned as a perturbation (RHOOUT) or as the full density ! (RHOFULL). Note that only the polynomial EOS choice will return ! a perturbation density; in other cases the full density is returned ! regardless of which argument is requested. ! ! This routine also computes derivatives of density with respect ! to temperature and salinity at level k from equation of state ! if requested (ie the optional arguments are present). ! ! If $k = kk$ are equal the density for level k is returned. ! If $k \neq kk$ the density returned is that for a parcel ! adiabatically displaced from level k to level kk. ! implicit none ! !INPUT PARAMETERS: integer, intent(in) :: & k, &! depth level index kk ! level to which water is adiabatically ! displaced real*8, dimension(i0,j0), intent(in) :: & tempk, &! temperature at level k saltk ! salinity at level k ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0), optional, intent(out) :: & rhoout, &! perturbation density of water rhofull, &! full density of water drhodt, &! derivative of density with respect to temperature drhods ! derivative of density with respect to salinity !EOP !BOC !----------------------------------------------------------------------- ! local variables: !----------------------------------------------------------------------- integer :: dep real*8, dimension(i0,j0) :: & TQ,SQ, &! adjusted T,S sqr,denomk, &! work arrays work1,work2,work3,work4 real*8 :: pp, pp2 ! temporary pressure scalars real*8 :: & mwjfnums0t0, & mwjfnums0t1, & mwjfnums0t2, & mwjfnums0t3, & mwjfnums1t0, & mwjfnums1t1, & mwjfnums2t0 real*8 :: & mwjfdens0t0, & mwjfdens0t1, & mwjfdens0t2, & mwjfdens0t3, & mwjfdens0t4, & mwjfdens1t0, & mwjfdens1t1, & mwjfdens1t3, & mwjfdensqt0, & mwjfdensqt2 !----------------------------------------------------------------------- ! valid ranges and pressure as function of depth !----------------------------------------------------------------------- real*8, dimension(k1) :: & tmin, tmax, &! valid temperature range for level k smin, smax, &! valid salinity range for level k pressz ! ref pressure (bars) at each level !----------------------------------------------------------------------- ! MWJF EOS coefficients !----------------------------------------------------------------------- !*** these constants will be used to construct the numerator !*** factor unit change (kg/m^3 -> g/cm^3) into numerator terms real*8, parameter :: p001 = 0.001 real*8, parameter :: & mwjfnp0s0t0 = 9.99843699e+2 * p001, & mwjfnp0s0t1 = 7.35212840e+0 * p001, & mwjfnp0s0t2 = -5.45928211e-2 * p001, & mwjfnp0s0t3 = 3.98476704e-4 * p001, & mwjfnp0s1t0 = 2.96938239e+0 * p001, & mwjfnp0s1t1 = -7.23268813e-3 * p001, & mwjfnp0s2t0 = 2.12382341e-3 * p001, & mwjfnp1s0t0 = 1.04004591e-2 * p001, & mwjfnp1s0t2 = 1.03970529e-7 * p001, & mwjfnp1s1t0 = 5.18761880e-6 * p001, & mwjfnp2s0t0 = -3.24041825e-8 * p001, & mwjfnp2s0t2 = -1.23869360e-11* p001 !*** these constants will be used to construct the denominator real*8, parameter :: & mwjfdp0s0t0 = 1.0e+0, & mwjfdp0s0t1 = 7.28606739e-3, & mwjfdp0s0t2 = -4.60835542e-5, & mwjfdp0s0t3 = 3.68390573e-7, & mwjfdp0s0t4 = 1.80809186e-10, & mwjfdp0s1t0 = 2.14691708e-3, & mwjfdp0s1t1 = -9.27062484e-6, & mwjfdp0s1t3 = -1.78343643e-10, & mwjfdp0sqt0 = 4.76534122e-6, & mwjfdp0sqt2 = 1.63410736e-9, & mwjfdp1s0t0 = 5.30848875e-6, & mwjfdp2s0t3 = -3.03175128e-16, & mwjfdp3s0t1 = -1.27934137e-17 !*** prevent problems with garbage on land points or ghost cells TQ = min(tempk, 1000.0_r8) TQ = max(TQ, -1000.0_r8) SQ = min(saltk, 1000.0_r8) SQ = max(saltk, 0.0_r8) !------------------------------------------------------------ do dep=1,k1 pressz(dep) = pressure(z(2*dep)*.01_r8) ! cm to m enddo pp = 10.0_r8*pressz(kk) SQ = 1000.0_r8*SQ !#ifdef CCSMCOUPLED ! call shr_vmath_sqrt(SQ, SQR, nx_block*ny_block) !#else SQR = sqrt(SQ) !#endif !*** !*** first calculate numerator of MWJF density [P_1(S,T,p)] !*** mwjfnums0t0 = mwjfnp0s0t0 + pp*(mwjfnp1s0t0 + pp*mwjfnp2s0t0) mwjfnums0t1 = mwjfnp0s0t1 mwjfnums0t2 = mwjfnp0s0t2 + pp*(mwjfnp1s0t2 + pp*mwjfnp2s0t2) mwjfnums0t3 = mwjfnp0s0t3 mwjfnums1t0 = mwjfnp0s1t0 + pp*mwjfnp1s1t0 mwjfnums1t1 = mwjfnp0s1t1 mwjfnums2t0 = mwjfnp0s2t0 WORK1 = mwjfnums0t0 + TQ * (mwjfnums0t1 + TQ * (mwjfnums0t2 + & mwjfnums0t3 * TQ)) + SQ * (mwjfnums1t0 + & mwjfnums1t1 * TQ + mwjfnums2t0 * SQ) !*** !*** now calculate denominator of MWJF density [P_2(S,T,p)] !*** mwjfdens0t0 = mwjfdp0s0t0 + pp*mwjfdp1s0t0 mwjfdens0t1 = mwjfdp0s0t1 + pp**3 * mwjfdp3s0t1 mwjfdens0t2 = mwjfdp0s0t2 mwjfdens0t3 = mwjfdp0s0t3 + pp**2 * mwjfdp2s0t3 mwjfdens0t4 = mwjfdp0s0t4 mwjfdens1t0 = mwjfdp0s1t0 mwjfdens1t1 = mwjfdp0s1t1 mwjfdens1t3 = mwjfdp0s1t3 mwjfdensqt0 = mwjfdp0sqt0 mwjfdensqt2 = mwjfdp0sqt2 WORK2 = mwjfdens0t0 + TQ * (mwjfdens0t1 + TQ * (mwjfdens0t2 + & TQ * (mwjfdens0t3 + mwjfdens0t4 * TQ))) + & SQ * (mwjfdens1t0 + TQ * (mwjfdens1t1 + TQ*TQ*mwjfdens1t3)+ & SQR * (mwjfdensqt0 + TQ*TQ*mwjfdensqt2)) DENOMK = 1.0_r8/WORK2 if (present(RHOOUT)) then RHOOUT = WORK1*DENOMK endif if (present(RHOFULL)) then RHOFULL = WORK1*DENOMK endif if (present(DRHODT)) then WORK3 = &! dP_1/dT mwjfnums0t1 + TQ * (2.0_r8*mwjfnums0t2 + & 3.0_r8*mwjfnums0t3 * TQ) + mwjfnums1t1 * SQ WORK4 = &! dP_2/dT mwjfdens0t1 + SQ * mwjfdens1t1 + & TQ * (2.0_r8*(mwjfdens0t2 + SQ*SQR*mwjfdensqt2) + & TQ * (3.0_r8*(mwjfdens0t3 + SQ * mwjfdens1t3) + & TQ * 4.0_r8*mwjfdens0t4)) DRHODT = (WORK3 - WORK1*DENOMK*WORK4)*DENOMK endif if (present(DRHODS)) then WORK3 = &! dP_1/dS mwjfnums1t0 + mwjfnums1t1 * TQ + 2.0_r8*mwjfnums2t0 * SQ WORK4 = mwjfdens1t0 + &! dP_2/dS TQ * (mwjfdens1t1 + TQ*TQ*mwjfdens1t3) + & 1.5_r8*SQR*(mwjfdensqt0 + TQ*TQ*mwjfdensqt2) DRHODS = (WORK3 - WORK1*DENOMK*WORK4)*DENOMK * 1000.0_r8 endif end subroutine state function pressure(depth) ! !INPUT PARAMETERS: real*8, intent(in) :: depth ! depth in meters ! !OUTPUT PARAMETERS: real*8 :: pressure ! pressure [bars] pressure = 0.059808_r8*(exp(-0.025_r8*depth) - 1.0_r8) & + 0.100766_r8*depth + 2.28405e-7_r8*depth**2 end function pressure subroutine ri_iwmix(dbloc, visc, vdc, uuu, vvv) ! !DESCRIPTION: ! Computes viscosity and diffusivity coefficients for the interior ! ocean due to shear instability (richardson number dependent), ! internal wave activity, and to static instability (Ri < 0). implicit none ! !INPUT PARAMETERS: real*8, dimension(i0,j0,k1), intent(in) :: & uuu ! U velocities at current time real*8, dimension(i0,j0,k1), intent(in) :: & vvv, &! V velocities at current time dbloc ! buoyancy difference between adjacent levels ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0,0:k0,2), intent(inout) :: & vdc ! diffusivity for tracer diffusion ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0,0:k0), intent(out) :: & visc ! viscosity !EOP !BOC !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & k, &! index for vertical levels i,j, &! horizontal loop indices n ! vertical smoothing index real*8, dimension(i0,j0) :: & vshear, &! (local velocity shear)^2 ri_loc, &! local Richardson number fri, &! function of Ri for shear work1 real*8, dimension(i0,j0) :: kvmix,kvmix_m !----------------------------------------------------------------------- ! compute mixing at each level !----------------------------------------------------------------------- kvmix = 0.0_r8 kvmix_m = 0.0_r8 visc(:,:,0) = 0.0_r8 do k = 1,k1 !----------------------------------------------------------------------- ! compute velocity shear squared and average to T points: ! VSHEAR = (UUU(k)-UUU(k+1))**2+(VVV(k)-VVV(k+1))**2 ! Use FRI here as a temporary. !----------------------------------------------------------------------- if (k < k1) then fri = (uuu(:,:,k)-uuu(:,:,k+1))**2 + & (vvv(:,:,k)-vvv(:,:,k+1))**2 vshear=fri else vshear = 0.0_r8 endif !----------------------------------------------------------------------- ! compute local richardson number ! use visc array as temporary Ri storage to be smoothed !----------------------------------------------------------------------- ri_loc = dbloc(:,:,k)*(zgrid(k)-zgrid(k+1))/(vshear + eps) visc(:,:,k) = merge( ri_loc, visc(:,:,k-1), k <= int(kb(:,:)) ) enddo !----------------------------------------------------------------------- ! vertically smooth Ri num_v_smooth_Ri times with 1-2-1 weighting ! result again stored temporarily in VISC and use RI_LOC and FRI ! as temps !----------------------------------------------------------------------- do n = 1,num_v_smooth_Ri fri = 0.250_r8 * visc(:,:,1) visc(:,:,k1+1) = visc(:,:,k1) do k=1,k1 do j=1,j0 do i=1,i0 ri_loc(i,j) = visc(i,j,k) if ( int(kb(i,j)) >= 1 ) then visc(i,j,k) = fri(i,j) + 0.500_r8*ri_loc(i,j) & + 0.250_r8*visc(i,j,k+1) endif fri(i,j) = 0.250_r8*ri_loc(i,j) enddo enddo enddo enddo !----------------------------------------------------------------------- ! now that we have a smoothed Ri field, finish computing coeffs ! at each level !----------------------------------------------------------------------- if ( ltidal_mixing ) tidal_diff(:,:,:) = 0.0_r8 ! write(*,*) 'kilmer_ri3' do k = 1,k1 !----------------------------------------------------------------------- ! ! if Ri-number mixing requested, ! evaluate function of Ri for shear instability: ! for 0 < Ri < Riinfty, function = (1 - (Ri/Riinfty)**2)**3 ! for Ri > Riinfty, function = 0 ! for Ri < 0 , function = 1 ! compute contribution due to shear instability ! VISC holds smoothed Ri at k, but replaced by real VISC ! ! otherwise only use iw ! convection is added later ! !----------------------------------------------------------------------- if ( ltidal_mixing ) then !----------------------------------------------------------------------- ! ! consider the internal wave mixing first. rich_mix is used as the ! upper limit for internal wave mixing coefficient. bckgrnd_vvc ! was already multiplied by Prandtl. ! ! NOTE: no partial_bottom_cell implementation at this time ! !----------------------------------------------------------------------- work1 = dbloc(:,:,k)/(zgrid(k) - zgrid(k+1)) where (work1 > 1.0_r8) tidal_diff(:,:,k) = tidal_coef(:,:,k)/work1 endwhere ! Notes: ! (1) this step breaks backwards compatibility ! (2) check for k>2 was added to if statement to avoid ! out of bounds access work1 = Prandtl*min(bckgrnd_vvc(:,:,k)/Prandtl & + tidal_diff(:,:,k), tidal_mix_max) if ( k < k1 ) then kvmix_m(:,:) = work1(:,:) endif if ( k < k1 ) then vdc(:,:,k,2) = min(bckgrnd_vdc(:,:,k) + tidal_diff(:,:,k), & tidal_mix_max) kvmix(:,:) = vdc(:,:,k,2) endif if (lrich) then fri = min((max(visc(:,:,k),0.0_r8))/Riinfty, 1.0_r8) visc(:,:,k) = work1 + rich_mix*(1.0_r8 - fri*fri)**3 if ( k < k1 ) then vdc(:,:,k,2) = vdc(:,:,k,2) + rich_mix*(1.0_r8 - fri*fri)**3 vdc(:,:,k,1) = vdc(:,:,k,2) endif else visc(:,:,k) = work1 if ( k < k1 ) then vdc(:,:,k,1) = vdc(:,:,k,2) endif endif else ! .not. ltidal_mixing ! write(*,*) 'kilmer_km1',k if ( k < k1 ) then kvmix(:,:) = bckgrnd_vdc(:,:,k) kvmix_m(:,:) = bckgrnd_vvc(:,:,k) endif if (lrich) then fri = min((max(visc(:,:,k),0.0_r8))/Riinfty, 1.0_r8) visc(:,:,k ) = bckgrnd_vvc(:,:,k) + & rich_mix*(1.0_r8 - fri*fri)**3 if ( k < k1 ) then vdc(:,:,k,2) = bckgrnd_vdc(:,:,k) + & rich_mix*(1.0_r8 - fri*fri)**3 vdc(:,:,k,1) = vdc(:,:,k,2) endif else visc(:,:,k ) = bckgrnd_vvc(:,:,k) if ( k < k1 ) then vdc(:,:,k,2) = bckgrnd_vdc(:,:,k) vdc(:,:,k,1) = vdc(:,:,k,2) endif endif endif ! ltidal_mixing !----------------------------------------------------------------------- ! set seafloor values to zero !----------------------------------------------------------------------- do j=1,j0 do i=1,i0 if ( k >= int(kb(i,j)) ) then visc(i,j,k ) = 0.0_r8 vdc (i,j,k,1) = 0.0_r8 vdc (i,j,k,2) = 0.0_r8 endif enddo enddo enddo !----------------------------------------------------------------------- ! fill extra coefficients for blmix !----------------------------------------------------------------------- visc(:,:,0 ) = 0.0_r8 vdc (:,:,0,:) = 0.0_r8 visc(:,:,k1+1 ) = 0.0_r8 vdc (:,:,k1+1,:) = 0.0_r8 end subroutine ri_iwmix subroutine init_vmix_kpp implicit none !----------------------------------------------------------------------- ! local coefficients !----------------------------------------------------------------------- integer :: i,j,k real*8 :: bckgrnd_vdc_psis ! PSI diffusivity in northern hemisphere real*8 :: bckgrnd_vdc_psin ! PSI diffusivity in southern hemisphere real*8, dimension(k1) :: zw ! vert dist from sfc to bottom of layer real*8, parameter :: bckgrnd_vdc1 = 0.1_r8 ! background diffusivity (Ledwell) real*8, parameter :: bckgrnd_vdc2 = 0.0 ! variation in diffusivity real*8, parameter :: bckgrnd_vdc_eq = 0.01 ! equatorial diffusivity (Gregg) real*8, parameter :: bckgrnd_vdc_psim = 0.13 ! Max. PSI induced diffusivity (MacKinnon) real*8, parameter :: bckgrnd_vdc_ban = 1.0 ! Banda Sea diffusivity (Gordon) real*8, parameter :: bckgrnd_vdc_dpth = 2500.0e02 ! depth at which diff equals vdc1 real*8, parameter :: bckgrnd_vdc_linv = 4.5e-05 ! inverse length for transition region ! tidal mixing real*8, dimension(i0,j0) :: & vertical_func, & work1,HT, & tidal_energy_flux real*8 :: coef,rho_fw real*8, parameter :: & local_mixing_fraction = 0.33, & mixing_efficiency = 0.2, & vertical_decay_scale = 500.0e02, & tidal_mix_max = 100.0 if (lhoriz_varying_bckgrnd) then k = 1 do j=1,j0 do i=1,i0 bckgrnd_vdc_psis= bckgrnd_vdc_psim*exp(-(0.4_r8*(ydeg(j)+28.9_r8))**2.) bckgrnd_vdc_psin= bckgrnd_vdc_psim*exp(-(0.4_r8*(ydeg(j)-28.9_r8))**2.) bckgrnd_vdc(i,j,k)=bckgrnd_vdc_eq+bckgrnd_vdc_psin+bckgrnd_vdc_psis if ( ydeg(j) .lt. -10.0_r8 ) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc(i,j,k) + bckgrnd_vdc1 elseif ( ydeg(j) .le. 10.0_r8 ) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc(i,j,k) + & bckgrnd_vdc1 * (ydeg(j)/10.0_r8)**2. else bckgrnd_vdc(i,j,k)=bckgrnd_vdc(i,j,k) + bckgrnd_vdc1 endif !---------------- ! North Banda Sea !---------------- if ( (ydeg(j) .lt. -1.0_r8) .and. (ydeg(j) .gt. -4.0_r8) .and. & (xdeg(i) .gt. 103.0_r8) .and. (xdeg(i) .lt. 134.0_r8)) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc_ban endif !----------------- ! Middle Banda Sea !----------------- if ( (ydeg(j) .le. -4.0_r8) .and. (ydeg(j) .gt. -7.0_r8) .and. & (xdeg(i) .gt. 106.0_r8) .and. (xdeg(i) .lt. 140.0_r8)) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc_ban endif !---------------- ! South Banda Sea !---------------- if ( (ydeg(j) .le. -7.0_r8) .and. (ydeg(j) .gt. -8.3_r8) .and. & (xdeg(i) .gt. 111.0_r8) .and. (xdeg(i) .lt. 142.0_r8)) then bckgrnd_vdc(i,j,k) = bckgrnd_vdc_ban endif bckgrnd_vvc(i,j,k) = Prandtl*bckgrnd_vdc(i,j,k) end do ! i end do ! j do k=2,k1 bckgrnd_vdc(:,:,k) = bckgrnd_vdc(:,:,1) bckgrnd_vvc(:,:,k) = bckgrnd_vvc(:,:,1) enddo else ! .not. lhoriz_varying_bckgrnd do k=1,k1 zw(k) = z(2*k+1) bckgrnd_vdc(:,:,k) = bckgrnd_vdc1 + bckgrnd_vdc2* & atan(bckgrnd_vdc_linv* & (zw(k)-bckgrnd_vdc_dpth)) bckgrnd_vvc(:,:,k) = Prandtl*bckgrnd_vdc(:,:,k) if (bckgrnd_vdc2 /= 0.0_r8 .and. myid == 0) then write (*,*)'bckgrnd_vdc=', bckgrnd_vdc(1,1,k) endif end do endif ! lhoriz_varying_bckgrnd ! initialize grid info zgrid(0) = eps hwide(0) = eps zgrid(k0) = -z(k1+k0) hwide(k0) = eps do k =1,k1 zgrid(k) = -z(2*k) hwide(k) = 1.0_r8/odz(k) enddo ! initialize cg and Vtc Vtc = sqrt(0.2_r8/c_s/epssfc)/vonkar**2 cg = cstar*vonkar*(c_s*vonkar*epssfc)**0.33_r8 ! stokes ratio do j = 1,j0 fstokes(:,j) = 11.0_r8 - max(5.0_r8*cos(3.0_r8*ydeg(j)) ,0.0_r8 ) enddo ! coriolis coeff fcort(1) = f(1) fcort(j0) = f(j2) do j=1,j2 fcort(j+1) = f(j) enddo if (ltidal_mixing) then tidal_diff = 0.0_r8 tidal_coef = 0.0_r8 rho_fw = 1.0_r8 HT = 0.0 do k = 1,k1 do j = 1,j0 do i = i,i0 if (k == int(kb(i,j))) HT(i,j) = z(2*k+1) enddo enddo enddo call read_tide(tidal_energy_file,tidal_energy_flux) write (*,*) ' tidal energy file read: ', trim(tidal_energy_file) !----------------------------------------------------------------------- ! convert TIDAL_ENERGY_FLUX from W/m^2 to gr/s^3. !----------------------------------------------------------------------- tidal_energy_flux(:,:) = 1000.0_r8 * tidal_energy_flux(:,:) !----------------------------------------------------------------------- ! compute the time independent part of the tidal mixing coefficients !----------------------------------------------------------------------- work1(:,:) = 0.0 do k=1,k1 where ( k < int(kb(:,:)) ) work1(:,:) = work1(:,:) + & exp(-(HT(:,:) - z(2*k+1))/vertical_decay_scale) / odzw(k+1) !dzw(k) endwhere end do ! k tidal_coef(:,:,:) = 0.0 coef = local_mixing_fraction * mixing_efficiency / rho_fw do k=1,k1 where ( k <= int(kb(:,:)) ) vertical_func(:,:) = & exp(-(HT(:,:) - z(2*k+1)) & /vertical_decay_scale) / work1(:,:) tidal_coef(:,:,k) = coef * & tidal_energy_flux(:,:) * vertical_func(:,:) elsewhere tidal_coef(:,:,k) = 0.0 endwhere enddo ! k end if ! ltidal_mixong end subroutine init_vmix_kpp subroutine ddmix(VDC, TRCR) ! !DESCRIPTION: ! $R_\rho$ dependent interior flux parameterization. ! Add double-diffusion diffusivities to Ri-mix values at blending ! interface and below. implicit none ! !INPUT PARAMETERS: real*8, dimension(i0,j0,k1,2), intent(in) :: & TRCR ! tracers at current time ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0,0:k0,2),intent(inout) :: & VDC ! diffusivity for tracer diffusion !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: k,kup,knxt real*8, dimension(i0,j0) :: & ALPHADT, &! alpha*DT across interfaces BETADS, &! beta *DS across interfaces RRHO, &! dd density ratio DIFFDD, &! dd diffusivity scale prtl_tmp ! prandtl number real*8, dimension(i0,j0,2) :: & TALPHA, &! temperature expansion coeff SBETA ! salinity expansion coeff !----------------------------------------------------------------------- ! compute alpha*DT and beta*DS at interfaces. use RRHO and ! PRANDTL for temporary storage for call to state !----------------------------------------------------------------------- kup = 1 knxt = 2 prtl_tmp = merge(-2.0_r8,TRCR(:,:,1,1),TRCR(:,:,1,1) < -2.0_r8) call state(1, 1, prtl_tmp, TRCR(:,:,1,2),RHOFULL=RRHO, & DRHODT=TALPHA(:,:,kup), DRHODS=SBETA(:,:,kup)) do k=1,k1 if ( k < k1 ) then prtl_tmp = merge(-2.0_r8,TRCR(:,:,k+1,1),TRCR(:,:,k+1,1) < -2.0_r8) call state(k+1, k+1, prtl_tmp, TRCR(:,:,k+1,2), & RHOFULL=RRHO, DRHODT=TALPHA(:,:,knxt), & DRHODS= SBETA(:,:,knxt)) ALPHADT = -0.5_r8*(TALPHA(:,:,kup) + TALPHA(:,:,knxt)) & *(TRCR(:,:,k,1) - TRCR(:,:,k+1,1)) BETADS = 0.5_r8*( SBETA(:,:,kup) + SBETA(:,:,knxt)) & *(TRCR(:,:,k,2) - TRCR(:,:,k+1,2)) kup = knxt knxt = 3 - kup else ALPHADT = 0.0_r8 BETADS = 0.0_r8 endif !----------------------------------------------------------------------- ! salt fingering case !----------------------------------------------------------------------- where ( ALPHADT > BETADS .and. BETADS > 0.0_r8 ) RRHO = MIN(ALPHADT/BETADS, Rrho0) DIFFDD = dsfmax*(1.0_r8-(RRHO-1.0_r8)/(Rrho0-1.0_r8))**3 VDC(:,:,k,1) = VDC(:,:,k,1) + 0.7_r8*DIFFDD VDC(:,:,k,2) = VDC(:,:,k,2) + DIFFDD endwhere !----------------------------------------------------------------------- ! diffusive convection !----------------------------------------------------------------------- where ( ALPHADT < 0.0_r8 .and. BETADS < 0.0_r8 .and. ALPHADT > BETADS ) RRHO = ALPHADT / BETADS DIFFDD = 1.5e-2_r8*0.909_r8* & exp(4.6_r8*exp(-0.54_r8*(1.0_r8/RRHO-1.0_r8))) prtl_tmp = 0.15_r8*RRHO elsewhere RRHO = 0.0_r8 DIFFDD = 0.0_r8 prtl_tmp = 0.0_r8 endwhere where (RRHO > 0.5_r8) prtl_tmp = (1.85_r8 - 0.85_r8/RRHO)*RRHO VDC(:,:,k,1) = VDC(:,:,k,1) + DIFFDD VDC(:,:,k,2) = VDC(:,:,k,2) + prtl_tmp*DIFFDD end do end subroutine ddmix subroutine bldepth(trans,dbloc, dbsfc, trcr, uuu, vvv, stf, shf_qsw,& hblt, ustar, bfsfc, stable, kbl, smf, smft) ! !DESCRIPTION: ! This routine computes the ocean boundary layer depth defined as ! the shallowest depth where the bulk Richardson number is equal to ! a critical value, Ricr. ! ! NOTE: bulk richardson numbers are evaluated by computing ! differences between values at zgrid(kl) $< 0$ and surface ! reference values. currently, the reference values are equal ! to the values in the surface layer. when using higher ! vertical grid resolution, these reference values should be ! computed as the vertical averages from the surface down to ! epssfc*zgrid(kl). ! ! This routine also computes where surface forcing is stable ! or unstable (STABLE) implicit none ! !INPUT PARAMETERS: real*8, dimension(i0,j0,k1,2), intent(in) :: & trcr ! tracers at current time real*8, dimension(i0,j0,k1), intent(in) :: & uuu,vvv, &! velocities at current time dbloc, &! buoyancy difference between adjacent levels dbsfc, &! buoyancy difference between level and surface trans ! transmission (1 to 0) from surface (calculated in kppglo) real*8, dimension(i0,j0,2), intent(in) :: & stf ! surface forcing for all tracers real*8, dimension(i0,j0), intent(in) :: & shf_qsw ! short-wave forcing real*8, dimension(i0,j0,2), intent(in), optional :: & smf, &! surface momentum forcing at U points smft ! surface momentum forcing at T points ! *** either one or the other (not ! *** both) should be passed ! !OUTPUT PARAMETERS: integer, dimension(i0,j0), intent(out) :: & kbl ! index of first lvl below hbl real*8, dimension(i0,j0), intent(out) :: & hblt, &! boundary layer depth bfsfc, &! Bo+radiation absorbed to d stable, &! =1 stable forcing; =0 unstab ustar ! surface friction velocity !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & i,j,k, &! loop indices kupper, kup, kdn, ktmp, kl ! vertical level indices real*8, dimension(i0,j0) :: & vshear, &! (velocity shear re sfc)^2 sigma, &! d/hbl wm, ws, &! turb vel scale functions bo, &! surface buoyancy forcing bosol, &! radiative buoyancy forcing talpha, &! temperature expansion coeff sbeta, &! salinity expansion coeff rho1, &! density at the surface work, &! temp array zkl, &! depth at current z level b_frqncy, &! buoyancy frequency rsh_hblt, &! resolved shear contribution to HBLT (fraction) hlangm, &! Langmuir depth hekman, &! Eckman depth limit hlimit ! limit to mixed-layer depth ! (= min(HEKMAN,HMONOB)) real*8, dimension(i0,j0,3) :: & ri_bulk, &! Bulk Ri number at 3 lvls hmonob ! Monin-Obukhov depth limit real*8 :: & absorb_frac, &! shortwave absorption frac sqrt_arg, &! dummy sqrt argument z_upper, z_up ! upper depths for RI_BULK interpolation real*8 :: & a_co, b_co, c_co ! coefficients of the quadratic equation ! $(a_{co}z^2+b_{co}|z|+c_{co}=Ri_b) used to ! find the boundary layer depth. when ! finding the roots, c_co = c_co - Ricr real*8 :: & slope_up ! slope of the above quadratic equation ! at zup. this is used as a boundary ! condition to determine the coefficients. !----------------------------------------------------------------------- ! compute friction velocity USTAR. !----------------------------------------------------------------------- if (present(smft)) then ustar = sqrt(sqrt(smft(:,:,1)**2 + smft(:,:,2)**2)) else work = sqrt(sqrt(smf(:,:,1)**2 + smf(:,:,2)**2)) ustar = work endif !----------------------------------------------------------------------- ! compute density and expansion coefficients at surface !----------------------------------------------------------------------- work = merge(-2.0_r8,trcr(:,:,1,1),trcr(:,:,1,1) < -2.0_r8) call state(1,1,work,trcr(:,:,1,2), & rhofull=rho1, drhodt=talpha, drhods=sbeta) !----------------------------------------------------------------------- ! compute turbulent and radiative sfc buoyancy forcing !----------------------------------------------------------------------- do j=1,j0 do i=1,i0 if (rho1(i,j) /= 0.0_r8) then bo (i,j) = grav*(-talpha(i,j)*stf(i,j,1) - & sbeta (i,j)*stf(i,j,2))/rho1(i,j) bosol(i,j) = -grav*talpha(i,j)*shf_qsw(i,j)/rho1(i,j) else bo (i,j) = 0.0_r8 bosol(i,j) = 0.0_r8 endif end do end do !----------------------------------------------------------------------- ! Find bulk Richardson number at every grid level until > Ricr ! max values when Ricr never satisfied are KBL = KMT and ! HBLT = -zgrid(KMT) ! ! NOTE: the reference depth is -epssfc/2.*zgrid(i,k), but the ! reference u,v,t,s values are simply the surface layer ! values and not the averaged values from 0 to 2*ref.depth, ! which is necessary for very fine grids(top layer < 2m ! thickness) ! ! ! Initialize hbl and kbl to bottomed out values ! Initialize HEKMAN and HLIMIT (= HMONOB until reset) to model bottom ! Initialize Monin Obukhov depth to value at z_up ! Set HMONOB=-zgrid(km) if unstable !----------------------------------------------------------------------- kupper = 1 kup = 2 kdn = 3 z_upper = 0.0_r8 z_up = zgrid(1) ri_bulk(:,:,kupper) = 0.0_r8 ri_bulk(:,:,kup) = 0.0_r8 kbl = merge( int(kb(:,:)), 1, (int(kb(:,:)) > 1) ) hlangm = 0.0_r8 do kl=1,k1 zkl = -zgrid(kl) do j=1,j0 do i=1,i0 if (kl == kbl(i,j)) hblt(i,j) = zkl(i,j) end do end do enddo if ( lcheckekmo ) then hekman = -zgrid(k1) + eps hlimit = -zgrid(k1) + eps if ( lshort_wave ) then bfsfc = bo + bosol*(1.0_r8-trans(:,:,1)) else bfsfc = bo endif stable = merge(1.0_r8, 0.0_r8, bfsfc >= 0.0_r8) bfsfc = bfsfc + stable*eps work = stable * cmonob*ustar*ustar*ustar/vonkar/bfsfc & + (stable -1.0_r8)*zgrid(k1) hmonob(:,:,kup) = merge( -z_up+eps, work, work <= -z_up ) endif rsh_hblt = 0.0_r8 !----------------------------------------------------------------------- ! compute velocity shear squared on U-grid and use the maximum ! of the four surrounding U-grid values for the T-grid. !----------------------------------------------------------------------- do kl = 2,k1 ! loop, 1846 work = (uuu(:,:,1)-uuu(:,:,kl))**2 + & (vvv(:,:,1)-vvv(:,:,kl))**2 ZKL = -zgrid(kl) VSHEAR(:,1) = 0.0_r8 do j=2,j0 vshear(1,j) = 0.0_r8 do i=2,i0 vshear(i,j) = max(work(i,j ), work(i-1,j ), & work(i,j-1), work(i-1,j-1)) enddo enddo !----------------------------------------------------------------------- ! compute bfsfc= Bo + radiative contribution down to hbf * hbl ! add epsilon to BFSFC to ensure never = 0 !----------------------------------------------------------------------- if (lshort_wave) then bfsfc = bo + bosol*(1.0_r8-trans(:,:,kl-1)) else bfsfc = bo endif stable = merge(1.0_r8, 0.0_r8, bfsfc >= 0.0_r8) bfsfc = bfsfc + stable*eps !----------------------------------------------------------------------- ! compute the Ekman and Monin Obukhov depths using above stability !----------------------------------------------------------------------- if (lcheckekmo) then fcort(1) = f(1) fcort(j0) = f(j2) do j=1,j2 fcort(j+1) = f(j) enddo do j=1,j0 do i=1,i0 if ( stable(i,j) > 0.5_r8 .and. hekman(i,j) >= -zgrid(k1) ) then hekman(i,j) = max(zkl(i,j), & cekman*ustar(i,j)/(abs(fcort(j))+eps)) endif end do end do hmonob(:,:,kdn) = stable*cmonob*ustar*ustar*ustar/vonkar/bfsfc + & (stable-1.0_r8)*zgrid(k1) do j=1,j0 do i=1,i0 if (hmonob(i,j,kdn) <= zkl(i,j) .and. & hmonob(i,j,kup) > -z_up) then work(i,j) = (hmonob(i,j,kdn) - hmonob(i,j,kup))/ & (z_up + zkl(i,j)) hlimit(i,j) = (hmonob(i,j,kdn) - work(i,j)*zkl(i,j))/ & (1.0_r8 - work(i,j)) endif end do end do endif !----------------------------------------------------------------------- ! compute velocity scales at sigma, for hbl = -zgrid(kl) !----------------------------------------------------------------------- sigma = epssfc call wscale(sigma, zkl, ustar, bfsfc, 2, wm, ws) !----------------------------------------------------------------------- ! compute the turbulent shear contribution to RI_BULK and store ! in WM. !----------------------------------------------------------------------- b_frqncy = sqrt(0.5_r8*(dbloc(:,:,kl) + abs(dbloc(:,:,kl)) + eps2)/ & (zgrid(kl)-zgrid(kl+1)) ) wm = zkl*ws*b_frqncy* & ( (Vtc/Ricr(kl))*max(2.1_r8 - 200.0_r8*b_frqncy,concv) ) !----------------------------------------------------------------------- ! compute bulk Richardson number at new level !----------------------------------------------------------------------- work = merge( (zgrid(1)-zgrid(kl))*dbsfc(:,:,kl), & 0.0_r8, int(kb(:,:)) >= kl) if ( linertial ) then ri_bulk(:,:,kdn) = work/(vshear+wm+ustar*bolus_sp(:,:)+eps) else ri_bulk(:,:,kdn) = work/(vshear+wm+eps) endif !----------------------------------------------------------------------- ! find hbl where Rib = Ricr. if possible, use a quadratic ! interpolation. if not, linearly interpolate. the quadratic ! equation coefficients are determined using the slope and ! Ri_bulk at z_up and Ri_bulk at zgrid(kl). the slope at ! z_up is computed linearly between z_upper and z_up. ! compute Langmuir depth always !----------------------------------------------------------------------- do j=1,j0 do i=1,i0 if ( kbl(i,j) == int(kb(i,j)) .and. & ri_bulk(i,j,kdn) > Ricr(kl) ) then slope_up = (ri_bulk(i,j,kupper) - ri_bulk(i,j,kup))/ & (z_up - z_upper) a_co = (ri_bulk(i,j,kdn) - ri_bulk(i,j,kup) - & slope_up*(zkl(i,j) + z_up) )/(z_up + zkl(i,j))**2 b_co = slope_up + 2.0_r8 * a_co * z_up c_co = ri_bulk(i,j,kup) + & z_up*(a_co*z_up + slope_up) - Ricr(kl) sqrt_arg = b_co**2 - 4.0_r8*a_co*c_co if ( ( abs(b_co) > eps .and. abs(a_co)/abs(b_co) <= eps ) & .or. sqrt_arg <= 0.0_r8 ) then hblt(i,j) = -z_up + (z_up + zkl(i,j)) * & (Ricr(kl) - ri_bulk(i,j,kup))/ & (ri_bulk(i,j,kdn) - ri_bulk(i,j,kup)) else hblt(i,j) = (-b_co + sqrt(sqrt_arg)) / (2.0_r8*a_co) endif kbl(i,j) = kl rsh_hblt(i,j) = (vshear(i,j)*Ricr(kl)/ & (dbsfc(i,j,kl)+eps))/hblt(i,j) hlangm(i,j) = ustar(i,j)*sqrt(fstokes(i,j)*zkl(i,j)/(dbsfc(i,j,kl)+eps) ) & / 0.9_r8 endif enddo enddo !----------------------------------------------------------------------- ! swap klevel indices and move to next level !----------------------------------------------------------------------- ktmp = kupper kupper = kup kup = kdn kdn = ktmp z_upper = z_up z_up = zgrid(kl) end do ! end loop, 1846 !----------------------------------------------------------------------- ! apply Langmuir parameterization if requested !----------------------------------------------------------------------- if ( llangmuir ) then do kl = k1,2,-1 where ( hlangm > hblt .and. & hlangm > -zgrid(kl-1) .and. & hlangm <= zkl ) hblt = hlangm kbl = kl end where enddo endif !----------------------------------------------------------------------- ! first combine Ekman and Monin-Obukhov depth limits. then apply ! these restrictions to HBLT. note that HLIMIT is set to -zgrid(km) ! in unstable forcing. !----------------------------------------------------------------------- if ( lcheckekmo ) then where ( hekman < hlimit ) hlimit = hekman do kl = 2,k1 where ( hlimit < hblt .and. & hlimit > -zgrid(kl-1) .and. & hlimit <= zkl ) hblt = hlimit kbl = kl end where enddo endif !----------------------------------------------------------------------- ! apply a Gaussian filter !----------------------------------------------------------------------- call smooth_hblt (.true., .false., hblt=hblt, kbl=kbl) !----------------------------------------------------------------------- ! correct stability and buoyancy forcing for SW up to boundary layer !----------------------------------------------------------------------- if (lshort_wave) bfsfc = bo + bosol*(1.0_r8-trans(:,:,1)) stable = merge(1.0_r8, 0.0_r8, bfsfc >= 0.0_r8) bfsfc = bfsfc + stable * eps ! ensures bfsfc never=0 end subroutine bldepth subroutine wscale(sigma, hbl, ustar, bfsfc, m_or_s, wm, ws) ! !DESCRIPTION: ! Computes turbulent velocity scales. ! ! For $\zeta \geq 0, ! w_m = w_s = \kappa U^\star/(1+5\zeta)$ ! ! For $\zeta_m \leq \zeta < 0, ! w_m = \kappa U^\star (1-16\zeta)^{1\over 4}$ ! ! For $\zeta_s \leq \zeta < 0, ! w_s = \kappa U^\star (1-16\zeta)^{1\over 2}$ ! ! For $\zeta < \zeta_m, ! w_m = \kappa U^\star (a_m - c_m\zeta)^{1\over 3}$ ! ! For $\zeta < \zeta_s, ! w_s = \kappa U^\star (a_s - c_s\zeta)^{1\over 3}$ ! ! where $\kappa$ is the von Karman constant. implicit none ! !INPUT PARAMETERS: integer, intent(in) :: & m_or_s ! flag =1 for wm only, 2 for ws, 3 for both real*8, dimension(i0,j0), intent(in) :: & sigma, &! normalized depth (d/hbl) hbl, &! boundary layer depth bfsfc, &! surface buoyancy forcing ustar ! surface friction velocity ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0), intent(out) :: & wm, &! turb velocity scales: momentum ws ! turb velocity scales: tracer !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: i,j ! dummy loop indices real*8, dimension(i0,j0) :: & zeta, &! d/L or sigma*hbl/L(monin-obk) zetah ! sigma*hbl*vonkar*BFSFC or ZETA = ZETAH/USTAR**3 !----------------------------------------------------------------------- ! compute zetah and zeta - surface layer is special case !----------------------------------------------------------------------- zetah = sigma*hbl*vonkar*bfsfc zeta = zetah/(ustar**3 + eps) !----------------------------------------------------------------------- ! compute velocity scales for momentum !----------------------------------------------------------------------- if (m_or_s == 1 .or. m_or_s == 3) then do j=1,j0 do i=1,i0 if (zeta(i,j) >= 0.0_r8) then ! stable region wm(i,j) = vonkar*ustar(i,j)/(1.0_r8 + 5.0_r8*zeta(i,j)) else if (zeta(i,j) >= zeta_m) then wm(i,j) = vonkar*ustar(i,j)*(1.0_r8 - 16.0_r8*zeta(i,j))**0.25_r8 else wm(i,j) = vonkar*(a_m*(ustar(i,j)**3)-c_m*zetah(i,j))**0.33_r8 endif end do end do endif !----------------------------------------------------------------------- ! compute velocity scales for tracers !----------------------------------------------------------------------- if (m_or_s == 2 .or. m_or_s == 3) then do j=1,j0 do i=1,i0 if (zeta(i,j) >= 0.0_r8) then ws(i,j) = vonkar*ustar(i,j)/(1.0_r8 + 5.0_r8*zeta(i,j)) else if (zeta(i,j) >= zeta_s) then ws(i,j) = vonkar*ustar(i,j)*sqrt(1.0_r8 - 16.0_r8*zeta(i,j)) else ws(i,j) = vonkar*(a_s*(ustar(i,j)**3)-c_s*zetah(i,j))**0.33_r8 endif end do end do endif end subroutine wscale subroutine smooth_hblt(overwrite_hblt, use_hmxl, & hblt, kbl, smooth_out) ! !DESCRIPTION: ! This subroutine uses a 1-1-4-1-1 Laplacian filter one time ! on HBLT or HMXL to reduce any horizontal two-grid-point noise. ! If HBLT is overwritten, KBL is adjusted after smoothing. implicit none ! !INPUT PARAMETERS: logical, intent(in) :: & overwrite_hblt, & ! if .true., HBLT is overwritten ! if .false., the result is returned in ! a dummy array use_hmxl ! if .true., smooth HMXL ! if .false., smooth HBLT ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0), optional, intent(inout) :: & hblt ! boundary layer depth integer, dimension(i0,j0), optional, intent(inout) :: & kbl ! index of first lvl below hbl ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0), optional, intent(out) :: & smooth_out ! optional output array containing the ! smoothened field if overwrite_hblt is false !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: & i, j, k real*8, dimension(i0,j0) :: & work1, work2 real*8 :: & ccc, ccw, cce, ccn, ccs, & ! averaging weights, * change to ccc and ccs ztmp ! temp for level depth !----------------------------------------------------------------------- ! perform one smoothing pass since we cannot do the necessary ! boundary updates for multiple passes. !----------------------------------------------------------------------- ! if ( use_hmxl ) then ! work2 = hmxl(:,:) ! else work2 = hblt ! endif work1 = work2 do j=2,j0-1 do i=2,i0-1 if ( int(kb(i,j)) /= 0 ) then ccw = 0.125_r8 cce = 0.125_r8 ccn = 0.125_r8 ccs = 0.125_r8 ccc = 0.5_r8 if ( int(kb(i-1,j)) == 0 ) then ccc = ccc + ccw ccw = 0.0_r8 endif if ( int(kb(i+1,j)) == 0 ) then ccc = ccc + cce cce = 0.0_r8 endif if ( int(kb(i,j-1)) == 0 ) then ccc = ccc + ccs ccs = 0.0_r8 endif if ( int(kb(i,j+1)) == 0 ) then ccc = ccc + ccn ccn = 0.0_r8 endif work2(i,j) = ccw * work1(i-1,j) & + cce * work1(i+1,j) & + ccs * work1(i,j-1) & + ccn * work1(i,j+1) & + ccc * work1(i,j) endif enddo enddo do k=1,k1 do j=2,j0-1 do i=2,i0-1 ztmp = -zgrid(k) if ( k == int(kb(i,j)) .and. work2(i,j) > ztmp ) then work2(i,j) = ztmp endif enddo enddo enddo if ( overwrite_hblt .and. .not.use_hmxl ) then hblt = work2 do k=1,k1 do j=2,j0-1 do i=2,i0-1 ztmp = -zgrid(k) if ( int(kb(i,j)) /= 0 .and. & ( hblt(i,j) > -zgrid(k-1) ) .and. & ( hblt(i,j) <= ztmp ) ) kbl(i,j) = k enddo enddo enddo else smooth_out = work2 endif end subroutine smooth_hblt subroutine blmix(visc, vdc, hblt, ustar, bfsfc, stable, & kbl, ghat) ! !DESCRIPTION: ! This routine computes mixing coefficients within boundary layer ! which depend on surface forcing and the magnitude and gradient ! of interior mixing below the boundary layer (matching). These ! quantities have been computed in other routines. ! ! Caution: if mixing bottoms out at hbl = -zgrid(km) then ! fictitious layer at km+1 is needed with small but finite width ! hwide(km+1). implicit none ! !INPUT/OUTPUT PARAMETERS: real*8, dimension(i0,j0,0:k0), intent(inout) :: & visc ! interior mixing coeff on input ! combined interior/bndy layer coeff output real*8, dimension(i0,j0,0:k0,2), intent(inout) :: & vdc ! diffusivity for tracer diffusion ! !INPUT PARAMETERS: integer, dimension(i0,j0), intent(in) :: & kbl ! index of first lvl below hbl real*8, dimension(i0,j0), intent(in) :: & hblt, & ! boundary layer depth bfsfc, & ! surface buoyancy forcing stable, & ! =1 stable forcing; =0 unstab ustar ! surface friction velocity ! !OUTPUT PARAMETERS: real*8, dimension(i0,j0,k1),intent(out) :: & ghat ! non-local mixing coefficient !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer :: & k,kp1, &! dummy k level index i,j ! horizontal indices integer, dimension(i0,j0) :: & kn ! klvl closest to HBLT real*8, dimension(i0,j0,k1,3) :: & blmc ! bndy layer mixing coefs real*8, dimension(i0,j0,3) :: & gat1, &! shape function at sigma=1 dat1, &! derivative of shape function dkm1 ! bndy layer difs at kbl-1 lvl real*8, dimension(i0,j0) :: & wm,ws, &! turbulent velocity scales casea, &! =1 in case A, =0 in case B sigma, &! normalized depth (d/hbl) visch, &! viscosity at hbl difth, &! temp diffusivity at hbl difsh, &! tracer diffusivity at hbl delhat, r, dvdzup, dvdzdn, & viscp, diftp, difsp, f1, & work1,work2 !----------------------------------------------------------------------- ! compute velocity scales at hbl !----------------------------------------------------------------------- sigma = epssfc call wscale(sigma, hblt, ustar, bfsfc, 3, wm, ws) !----------------------------------------------------------------------- ! determine caseA = 0 if closer to KBL than KBL-1 ! KN is then the closest klevel to HBLT !----------------------------------------------------------------------- do j=1,j0 do i=1,i0 k = kbl(i,j) casea(i,j) = 0.5_r8 + SIGN(0.5_r8, -zgrid(k)-0.5*hwide(k)-hblt(i,j)) enddo enddo kn = NINT(casea)*(kbl-1) + (1-NINT(casea))*kbl !----------------------------------------------------------------------- ! find the interior viscosities and derivatives at hbl by ! interpolating derivative values at vertical interfaces. compute ! matching conditions for shape function. !----------------------------------------------------------------------- f1 = stable*5.0_r8*bfsfc/(ustar**4+eps) do k=1,k1 do j=1,j0 do i=1,i0 if (k == kn(i,j)) then delhat(i,j) = 0.5_r8*hwide(k) - zgrid(k) - hblt(i,j) r (i,j) = 1.0_r8 - delhat(i,j) / hwide(k) dvdzup(i,j) = (visc(i,j,k-1) - visc(i,j,k ))/hwide(k) dvdzdn(i,j) = (visc(i,j,k ) - visc(i,j,k+1))/hwide(k+1) viscp (i,j) = 0.5_r8*( (1.0_r8-r(i,j))* & (dvdzup(i,j) + abs(dvdzup(i,j))) + & r(i,j) * & (dvdzdn(i,j) + abs(dvdzdn(i,j))) ) dvdzup(i,j) = (vdc(i,j,k-1,2) - vdc(i,j,k ,2))/hwide(k) dvdzdn(i,j) = (vdc(i,j,k ,2) - vdc(i,j,k+1,2))/hwide(k+1) difsp (i,j) = 0.5_r8*( (1.0_r8-r(i,j))* & (dvdzup(i,j) + abs(dvdzup(i,j))) + & r(i,j) * & (dvdzdn(i,j) + abs(dvdzdn(i,j))) ) dvdzup(i,j) = (vdc(i,j,k-1,1) - vdc(i,j,k ,1))/hwide(k) dvdzdn(i,j) = (vdc(i,j,k ,1) - vdc(i,j,k+1,1))/hwide(k+1) diftp (i,j) = 0.5_r8*( (1.0_r8-r(i,j))* & (dvdzup(i,j) + abs(dvdzup(i,j))) + & r(i,j) * & (dvdzdn(i,j) + abs(dvdzdn(i,j))) ) visch(i,j) = visc(i,j,k) + viscp(i,j)*delhat(i,j) difsh(i,j) = vdc(i,j,k,2) + difsp(i,j)*delhat(i,j) difsh(i,j) = vdc(i,j,k,1) + diftp(i,j)*delhat(i,j) gat1(i,j,1) = visch(i,j) / hblt(i,j) /(wm(i,j)+eps) dat1(i,j,1) = -viscp(i,j)/(wm(i,j)+eps) + & f1(i,j)*visch(i,j) gat1(i,j,2) = difsh(i,j) / hblt(i,j) /(ws(i,j)+eps) dat1(i,j,2) = -difsp(i,j)/(ws(i,j)+eps) + & f1(i,j)*difsh(i,j) gat1(i,j,3) = difth(i,j) / hblt(i,j) /(ws(i,j)+eps) dat1(i,j,3) = -diftp(i,j)/(ws(i,j)+eps) + & f1(i,j)*difth(i,j) endif end do end do enddo dat1 = min(dat1,0.0_r8) !----------------------------------------------------------------------- ! compute the dimensionless shape functions and diffusivities ! at the grid interfaces. also compute function for non-local ! transport term (GHAT). !----------------------------------------------------------------------- do k = 1,k1 sigma = (-zgrid(k) + 0.5_r8*hwide(k)) / hblt f1 = min(sigma,epssfc) call wscale(f1, hblt, ustar, bfsfc, 3, wm, ws) do j=1,j0 do i=1,i0 blmc(i,j,k,1) = hblt(i,j)*wm(i,j)*sigma(i,j)* & (1.0_r8 + sigma(i,j)*((sigma(i,j)-2.0_r8) + & (3.0_r8-2.0_r8*sigma(i,j))*gat1(i,j,1) + & (sigma(i,j)-1.0_r8)*dat1(i,j,1))) blmc(i,j,k,2) = hblt(i,j)*ws(i,j)*sigma(i,j)* & (1.0_r8 + sigma(i,j)*((sigma(i,j)-2.0_r8) + & (3.0_r8-2.0_r8*sigma(i,j))*gat1(i,j,2) + & (sigma(i,j)-1.0_r8)*dat1(i,j,2))) blmc(i,j,k,3) = hblt(i,j)*ws(i,j)*sigma(i,j)* & (1.0_r8 + sigma(i,j)*((sigma(i,j)-2.0_r8) + & (3.0_r8-2.0_r8*sigma(i,j))*gat1(i,j,3) + & (sigma(i,j)-1.0_r8)*dat1(i,j,3))) ghat(i,j,k) = (1.0_r8-stable(i,j))* cg/(ws(i,j)*hblt(i,j) +eps) end do end do end do !----------------------------------------------------------------------- ! find diffusivities at kbl-1 grid level !----------------------------------------------------------------------- do j=1,j0 do i=1,i0 k = kbl(i,j) - 1 sigma(i,j) = -zgrid(k)/hblt(i,j) enddo enddo f1 = min(sigma,epssfc) call wscale(f1, hblt, ustar, bfsfc, 3, wm, ws) do j=1,j0 do i=1,i0 dkm1(i,j,1) = hblt(i,j)*wm(i,j)*sigma(i,j)* & (1.0_r8+sigma(i,j)*((sigma(i,j)-2.0_r8) + & (3.0_r8-2.0_r8*sigma(i,j))*gat1(i,j,1) + & (sigma(i,j)-1.0_r8)*dat1(i,j,1))) dkm1(i,j,2) = hblt(i,j)*ws(i,j)*sigma(i,j)* & (1.0_r8+sigma(i,j)*((sigma(i,j)-2.0_r8) + & (3.0_r8-2.0_r8*sigma(i,j))*gat1(i,j,2) + & (sigma(i,j)-1.0_r8)*dat1(i,j,2))) dkm1(i,j,3) = hblt(i,j)*ws(i,j)*sigma(i,j)* & (1.0_r8+sigma(i,j)*((sigma(i,j)-2.0_r8) + & (3.0_r8-2.0_r8*sigma(i,j))*gat1(i,j,3) + & (sigma(i,j)-1.0_r8)*dat1(i,j,3))) end do end do !----------------------------------------------------------------------- ! compute the enhanced mixing !----------------------------------------------------------------------- do k=1,k1-1 where (k == (kbl - 1)) & delhat = (hblt + zgrid(k))/(zgrid(k)-zgrid(k+1)) do j=1,j0 do i=1,i0 if (k == (kbl(i,j) - 1)) then blmc(i,j,k,1) = (1.0_r8-delhat(i,j))*visc(i,j,k) + & delhat(i,j) *( & (1.0_r8-delhat(i,j))**2 *dkm1(i,j,1) + & delhat(i,j)**2 *(casea(i,j)*visc(i,j,k) + & (1.0_r8-casea(i,j))*blmc(i,j,k,1))) blmc(i,j,k,2) = (1.0_r8-delhat(i,j))*vdc(i,j,k,2) + & delhat(i,j)*( & (1.0_r8-delhat(i,j))**2 *dkm1(i,j,2) + & delhat(i,j)**2 *(casea(i,j)*vdc(i,j,k,2) + & (1.0_r8-casea(i,j))*blmc(i,j,k,2))) blmc(i,j,k,3) = (1.0_r8-delhat(i,j))*vdc(i,j,k,1) + & delhat(i,j) *( & (1.0_r8-delhat(i,j))**2 *dkm1(i,j,3) + & delhat(i,j)**2 *(casea(i,j)*vdc(i,j,k,1) + & (1.0_r8-casea(i,j))*blmc(i,j,k,3))) ghat(i,j,k) = (1.0_r8-casea(i,j)) * ghat(i,j,k) endif end do end do end do !----------------------------------------------------------------------- ! combine interior and boundary layer coefficients and nonlocal term !----------------------------------------------------------------------- do k=1,k1 do j=1,j0 do i=1,i0 if (k < kbl(i,j)) then visc(i,j,k) = blmc(i,j,k,1) vdc(i,j,k,2) = blmc(i,j,k,2) vdc(i,j,k,1) = blmc(i,j,k,3) else ghat(i,j,k) = 0.0_r8 endif end do end do enddo end subroutine blmix end module physics