module physics
    use const
    use comm
    use var 
    use io
implicit real*8 (a-h,o-z)
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 kppglo
end subroutine kppglo
!--------------------------------------------------------------------------


! ----------------------------------------------------------------------
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 qsurfglo(t,dt,odz)
    dimension t(i0,j0),odz(*)
    ! 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.
    ! cp for water is 2.5e4 ergs/gram/deg c.
    !check this cp value!!!
    tmp=dt*odz(1)/2.5e4
    dtemp=tmp*qdot
    do j=1,j2
	do i=1,i2
	    t(i+1,j+1)=t(i+1,j+1)+dtemp
	enddo
    enddo
end subroutine

! ----------------------------------------------------------------------
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


end module
